﻿// --------------------------------------------------------------------------------------------------------------------
// <copyright file="JDelone.cs" company="">
//   
// </copyright>
// <summary>
//   The face walk.
// </summary>
// --------------------------------------------------------------------------------------------------------------------

namespace FrontalMethod.Java
{
    using System;
    using System.Collections.Generic;
    using System.Diagnostics;
    using System.Drawing;
    using System.Drawing.Drawing2D;

    /**
 * A helper class for returning the result of a face walk.
 *
 * @author Mark Howison
 */

    /// <summary>
    /// The face walk.
    /// </summary>
    public class FaceWalk
    {
        /// <summary>
        /// The failed.
        /// </summary>
        public static int FAILED = -1;

        /// <summary>
        /// The coincident.
        /// </summary>
        public static int COINCIDENT;

        /// <summary>
        /// The clockwise.
        /// </summary>
        public static int CLOCKWISE = 1;

        /// <summary>
        /// The he.
        /// </summary>
        public HalfEdge he;

        /// <summary>
        /// The status.
        /// </summary>
        public int status;

        /// <summary>
        /// Initializes a new instance of the <see cref="FaceWalk"/> class.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        /// <param name="status">
        /// The status.
        /// </param>
        public FaceWalk(HalfEdge he, int status)
        {
            this.he = he;
            this.status = status;
        }
    }

    /**
 * Represents a boundary point, which has additional fields <tt>next</tt>
 * and <tt>prev</tt> to maintain a doubley-linked list of the boundary.
 * The <tt>side</tt> field can be used for encoding symmetry information
 * about the boundary.
 *
 * @author Mark Howison
 */

    /// <summary>
    /// The b point.
    /// </summary>
    public class BPoint : Point
    {
        /// <summary>
        /// The next.
        /// </summary>
        public BPoint next;

        /// <summary>
        /// The prev.
        /// </summary>
        public BPoint prev;

        /// <summary>
        /// The side.
        /// </summary>
        public int side;

        /**
         * Creates a point with 'boundary' type.
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="BPoint"/> class.
        /// </summary>
        public BPoint()
        {
            // super();
            this.type = BOUNDARY;
        }

        /**
         * Sets the 2D coordinates to those of <tt>p</tt>.
         *
         * @param p
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="BPoint"/> class.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        public BPoint(Point2d p)
            : base(p)
        {
            // super(p);
            this.type = BOUNDARY;
        }

        /**
         * Sets the 2D coordinates to <tt>(x,y)</tt>.
         * 
         * @param x
         * @param y
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="BPoint"/> class.
        /// </summary>
        /// <param name="x">
        /// The x.
        /// </param>
        /// <param name="y">
        /// The y.
        /// </param>
        public BPoint(double x, double y)
            : base(x, y)
        {
            // super(x, y);
            this.type = BOUNDARY;
        }

        /**
         * Creates a shallow copy of <tt>bp</tt>, copying the pointers and not
         * the objects.
         *
         * @param bp
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="BPoint"/> class.
        /// </summary>
        /// <param name="bp">
        /// The bp.
        /// </param>
        public BPoint(BPoint bp)
            : base(bp)
        {
            // super(bp);
            this.type = BOUNDARY;
            this.side = bp.side;
            this.pair = bp.pair;
            this.next = bp.next;
            this.prev = bp.prev;
        }

        /**
         * Up-casts this point's pair to a <tt>BPoint</tt>.
         *
         * @return
         */

        /// <summary>
        /// The get pair.
        /// </summary>
        /// <returns>
        /// </returns>
        public BPoint getPair()
        {
            return (BPoint)this.pair;
        }
    }

    /// <summary>
    /// 2D Point with double precision
    /// </summary>
    public class Point2d
    {
        /// <summary>
        /// X component of point
        /// </summary>
        protected double _X;

        /// <summary>
        /// Y component of point
        /// </summary>
        protected double _Y;

        /// <summary>
        /// Point number
        /// </summary>
        protected int number;

        /// <summary>
        /// Initializes a new instance of the <see cref="Point2d"/> class.
        /// </summary>
        public Point2d()
        {
            this._X = 0;
            this._Y = 0;
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="Point2d"/> class.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        public Point2d(Point2d p)
        {
            this._X = p.X;
            this._Y = p.Y;
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="Point2d"/> class. 
        /// Initializes a new instance of a point
        /// </summary>
        /// <param name="x">
        /// </param>
        /// <param name="y">
        /// </param>
        public Point2d(double x, double y)
        {
            this._X = x;
            this._Y = y;
        }

        /// <summary>
        /// Gets or sets the X component of the point
        /// </summary>
        public double X
        {
            get
            {
                return this._X;
            }

            set
            {
                this._X = value;
            }
        }

        /// <summary>
        /// Gets or sets the Y component of the point
        /// </summary>
        public double Y
        {
            get
            {
                return this._Y;
            }

            set
            {
                this._Y = value;
            }
        }

        /// <summary>
        /// Gets or sets Number.
        /// </summary>
        public int Number
        {
            get
            {
                return this.number;
            }

            set
            {
                this.number = value;
            }
        }

        /// <summary>
        /// Makes a planar checks for if the points is spatially equal to another point.
        /// </summary>
        /// <param name="other">
        /// Point to check against
        /// </param>
        /// <returns>
        /// True if X and Y values are the same
        /// </returns>
        public bool Equals2D(Point other)
        {
            return this.X == other.X && this.Y == other.Y;
        }
    }

    /**
     * Represents a 2D vertex in a triangulation. Contains references to a halfedge
     * in the star of the vertex, to 3D coordinates (for 3D meshing), and to
     * a "paired" point (for encoding symmetry information).
     *
     * @author Mark Howison
     */

    /// <summary>
    /// The point.
    /// </summary>
    public class Point : Point2d
    {
        /// <summary>
        /// The interior.
        /// </summary>
        public static int INTERIOR;

        /// <summary>
        /// The boundary.
        /// </summary>
        public static int BOUNDARY = 1;

        /// <summary>
        /// The deleted.
        /// </summary>
        public static int DELETED = 2;

        /// <summary>
        /// The translated.
        /// </summary>
        public static int TRANSLATED = 3;

        /// <summary>
        /// The xsection.
        /// </summary>
        public static int XSECTION = 4;

        /// <summary>
        /// The bounds.
        /// </summary>
        public static int BOUNDS = 5;

        /// <summary>
        /// The colors.
        /// </summary>
        public static Color[] COLORS = { Color.Magenta, Color.Cyan, Color.Red, Color.Orange, Color.Green, Color.Cyan };

        /// <summary>
        /// The he.
        /// </summary>
        public HalfEdge he;

        // public Coords3d coords3d = new Coords3d();
        /// <summary>
        /// The pair.
        /// </summary>
        public Point pair;

        /// <summary>
        /// The type.
        /// </summary>
        public int type = INTERIOR;

        /**
         * Constructs a new point with 2D coordinates (0,0).
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="Point"/> class.
        /// </summary>
        public Point()
        {
            // super();
        }

        /**
         * Constructs a new point with the same 2D coordinates as <tt>p</tt>.
         *
         * @param p
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="Point"/> class.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        public Point(Point2d p)
            : base(p)
        {
            // super(p);
        }

        /**
         * Constructs a new point with 2D coordinates <tt>(x,y)</tt>.
         *
         * @param p
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="Point"/> class.
        /// </summary>
        /// <param name="x">
        /// The x.
        /// </param>
        /// <param name="y">
        /// The y.
        /// </param>
        public Point(double x, double y)
            : base(x, y)
        {
            // super(x, y);
        }

        /**
         * Shallow copy of point <tt>p</tt>. All pointers are copied, but not
         * the underlying objects. 2D and 3D coordinates are set to those of
         * <tt>p</tt>.
         *
         * @param p
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="Point"/> class.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        public Point(Point p)
            : base(p)
        {
            // super(p);
            this.he = p.he;
            this.type = p.type;

            // coords3d.set(p.coords3d);
            this.pair = p.pair;
        }

        /**
         * Sets the point type (interior, boundary, etc.). Use the static int flags
         * to specify type.
         *
         * @param type
         */

        /// <summary>
        /// The set type.
        /// </summary>
        /// <param name="type">
        /// The type.
        /// </param>
        public void setType(int type)
        {
            this.type = type;
        }

        /**
         * Tests if this point is of <tt>type</tt>. Use the static int flags
         * to specify type.
         *
         * @param type
         * @return
         */

        /// <summary>
        /// The is type.
        /// </summary>
        /// <param name="type">
        /// The type.
        /// </param>
        /// <returns>
        /// The is type.
        /// </returns>
        public bool isType(int type)
        {
            if (this.type == type)
            {
                return true;
            }

            if (type == INTERIOR && this.type == TRANSLATED)
            {
                return true;
            }

            if (type == BOUNDARY && this.type == XSECTION)
            {
                return true;
            }

            return false;
        }

        /**
         * Gets the color associated with this point's type.
         *
         * @return
         */

        /// <summary>
        /// The get color.
        /// </summary>
        /// <returns>
        /// </returns>
        public Color getColor()
        {
            return COLORS[this.type];
        }
    }

    /**
 * A helper class that pairs the endpoints of an edge.
 *
 * @author Mark Howison
 */

    /// <summary>
    /// The edge.
    /// </summary>
    public class Edge
    {
        /// <summary>
        /// The p 1.
        /// </summary>
        public Point p1;

        /// <summary>
        /// The p 2.
        /// </summary>
        public Point p2;

        /// <summary>
        /// The type.
        /// </summary>
        public int type;

        /**
         * Constructs an edge with endpoints <tt>p1</tt> and <tt>p2</tt>.
         *
         * @param p1
         * @param p2
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="Edge"/> class.
        /// </summary>
        /// <param name="p1">
        /// The p 1.
        /// </param>
        /// <param name="p2">
        /// The p 2.
        /// </param>
        public Edge(Point p1, Point p2)
        {
            this.p1 = p1;
            this.p2 = p2;
        }

        /**
         * Copy constructor for an existing edge <tt>e</tt>. The pointers to the
         * endpoints are copied, but not the points themsleves (shallow copy).
         *
         * @param e
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="Edge"/> class.
        /// </summary>
        /// <param name="e">
        /// The e.
        /// </param>
        public Edge(Edge e)
        {
            this.p1 = e.p1;
            this.p2 = e.p2;
        }

        /**
         * Returns <tt>true</tt> if edges <tt>e1</tt> and <tt>e2</tt> have any
         * endpoint in common.
         *
         * @param e1
         * @param e2
         * @return
         */

        /// <summary>
        /// The adjacent.
        /// </summary>
        /// <param name="e1">
        /// The e 1.
        /// </param>
        /// <param name="e2">
        /// The e 2.
        /// </param>
        /// <returns>
        /// The adjacent.
        /// </returns>
        public static bool adjacent(Edge e1, Edge e2)
        {
            if (e1.p1 == e2.p1)
            {
                return true;
            }

            if (e1.p1 == e2.p2)
            {
                return true;
            }

            if (e1.p2 == e2.p1)
            {
                return true;
            }

            if (e1.p2 == e2.p2)
            {
                return true;
            }

            return false;
        }

        /**
         * Returns <tt>true</tt> if this edge has the same endpoints as <tt>e</tt>.
         *
         * @param e
         * @return
         */

        /// <summary>
        /// The equals.
        /// </summary>
        /// <param name="e">
        /// The e.
        /// </param>
        /// <returns>
        /// The equals.
        /// </returns>
        public bool equals(Edge e)
        {
            if (e.p1 == this.p1 && e.p2 == this.p2)
            {
                return true;
            }

            if (e.p1 == this.p2 && e.p2 == this.p1)
            {
                return true;
            }

            return false;
        }
    }

    /**
 * A halfedge is a directed segment that represents one side of an edge in
 * the triangulation. Each halfedge has pointers to an origin point,
 * the next halfedge, and its "sibling" halfedge, which points in the opposite
 * direction but shares the same edge of the triangulation.
 *
 * A valid triangulation will contain a 3-cycle of halfedges for each face,
 * every interior edge will correspond to a sibling pair of halfeldges, and
 * every boundary edge will correspond to a single halfedge with a null
 * sibling pointer.
 *
 * @author Mark Howison
 */

    /// <summary>
    /// The half edge.
    /// </summary>
    public class HalfEdge
    {
        /// <summary>
        /// The auxilary.
        /// </summary>
        public static int AUXILARY;

        /// <summary>
        /// The boundary.
        /// </summary>
        public static int BOUNDARY = 1;

        /// <summary>
        /// The constraint.
        /// </summary>
        public static int CONSTRAINT = 2;

        /// <summary>
        /// The flags.
        /// </summary>
        public static int FLAGS = 3;

        /// <summary>
        /// The fla g_ algorithm.
        /// </summary>
        public static int FLAG_ALGORITHM;

        /// <summary>
        /// The fla g_ draw.
        /// </summary>
        public static int FLAG_DRAW = 1;

        /// <summary>
        /// The fla g_ read.
        /// </summary>
        public static int FLAG_READ = 2;

        /// <summary>
        /// The colors.
        /// </summary>
        public static Color[] COLORS = { Color.Gray, Color.Blue, Color.Magenta };

        /// <summary>
        /// The flagged.
        /// </summary>
        private readonly bool[] flagged = new bool[FLAGS];

        /// <summary>
        /// The next.
        /// </summary>
        public HalfEdge next;

        /// <summary>
        /// The origin.
        /// </summary>
        public Point origin;

        /// <summary>
        /// The sibling.
        /// </summary>
        public HalfEdge sibling;

        /// <summary>
        /// The type.
        /// </summary>
        private int type = AUXILARY;

        /**
         * Constructs a halfedge with null pointers.
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="HalfEdge"/> class.
        /// </summary>
        public HalfEdge()
        {
        }

        /**
         * Constructs a halfedge with origin <tt>p</tt>, and null next and sibling
         * pointers.
         *
         * @param p
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="HalfEdge"/> class.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        public HalfEdge(Point p)
        {
            this.origin = p;
        }

        /**
         * Constructs a halfedge with origin <tt>p</tt>, type <tt>t</tt>, and
         * null next and sibling pointers.
         *
         * @param p
         * @param t
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="HalfEdge"/> class.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <param name="t">
        /// The t.
        /// </param>
        public HalfEdge(Point p, int t)
        {
            this.origin = p;
            this.type = t;
        }

        /**
         * Shallow copies <tt>he</tt>. All pointers are copied, but not the
         * underlying objects.
         *
         * @param he
         */

        /// <summary>
        /// Initializes a new instance of the <see cref="HalfEdge"/> class.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        public HalfEdge(HalfEdge he)
        {
            this.origin = he.origin;
            this.next = he.next;
            this.sibling = he.sibling;
            this.type = he.type;
            this.flagged = he.flagged;
        }

        /**
         * Gets the next halfedge.
         *
         * @return
         */

        /// <summary>
        /// The get next.
        /// </summary>
        /// <returns>
        /// </returns>
        public HalfEdge getNext()
        {
            return this.next;
        }

        /**
         * Gets the previous halfedge whose next pointer points to this halfedge.
         *
         * NOTE: In a valid triangulation, this could be accomplished by following
         * three next pointers. However, this method is robust for non-triangular
         * regions, and instead uses a potentially infinite while loop.
         *
         * @return
         */

        /// <summary>
        /// The get prev.
        /// </summary>
        /// <returns>
        /// </returns>
        public HalfEdge getPrev()
        {
            HalfEdge hePrev = this.next;
            while (true)
            {
                if (hePrev.next == this)
                {
                    break;
                }

                hePrev = hePrev.next;
            }

            return hePrev;
        }

        /**
         * Gets the color associated with this halfedge's type.
         *
         * @return
         */

        /// <summary>
        /// The get color.
        /// </summary>
        /// <returns>
        /// </returns>
        public Color getColor()
        {
            return COLORS[this.type];
        }

        /**
         * Gets this halfedge's type (interior, boundary, constrained etc.). Use
         * the static int flags to specify type.
         *
         * @return
         */

        /// <summary>
        /// The get type.
        /// </summary>
        /// <returns>
        /// The get type.
        /// </returns>
        public int getType()
        {
            return this.type;
        }

        /**
         * Sets the next halfedge pointer to <tt>he</tt>.
         *
         * @param he
         */

        /// <summary>
        /// The set next.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        public void setNext(HalfEdge he)
        {
            this.next = he;
            if (this.type != BOUNDARY)
            {
                if (this.next.origin != this.sibling.origin)
                {
                    ;
                }
            }
        }

        /**
         * Sets the sibling halfedge pointer to <tt>he</tt>.
         *
         * @param he
         */

        /// <summary>
        /// The set sibling.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        public void setSibling(HalfEdge he)
        {
            this.sibling = he;
        }

        /**
         * Sets this halfedge's type (interior, boundary, constrained etc.). Use
         * the static int flags to specify type.
         *
         * @param type
         */

        /// <summary>
        /// The set type.
        /// </summary>
        /// <param name="type">
        /// The type.
        /// </param>
        public void setType(int type)
        {
            this.type = type;
            if (this.sibling != null)
            {
                this.sibling.type = type;
            }
        }

        /**
         * Tests if this halfedges is of <tt>type</tt>. Use the static int flags
         * to specify type.
         *
         * @param type
         * @return
         */

        /// <summary>
        /// The is type.
        /// </summary>
        /// <param name="type">
        /// The type.
        /// </param>
        /// <returns>
        /// The is type.
        /// </returns>
        public bool isType(int type)
        {
            if (this.type == type)
            {
                return true;
            }
            else
            {
                return false;
            }
        }

        /**
         * Tests whether this halfedge's <tt>flag</tt> is set to true.
         *
         * Flags can be used to mark which edges or faces of a triangulation
         * have been already processed by a procedure or algorithm.
         *
         * @param flag
         * @return
         */

        /// <summary>
        /// The is flagged.
        /// </summary>
        /// <param name="flag">
        /// The flag.
        /// </param>
        /// <returns>
        /// The is flagged.
        /// </returns>
        public bool isFlagged(int flag)
        {
            return this.flagged[flag];
        }

        /**
         * Sets this halfedge's type to CONSTRAINT.
         */

        /// <summary>
        /// The constrain.
        /// </summary>
        public void constrain()
        {
            this.type = CONSTRAINT;
        }

        /**
         * Sets this halfedge's <tt>flag</tt> to true. Does not change the
         * sibling's flag.
         *
         * @param flag
         */

        /// <summary>
        /// The flag.
        /// </summary>
        /// <param name="flag">
        /// The flag.
        /// </param>
        public void flag(int flag)
        {
            this.flagged[flag] = true;
        }

        /**
         * Sets this halfedge and its sibling's <tt>flag</tt> to true.
         *
         * @param flag
         */

        /// <summary>
        /// The flag edge.
        /// </summary>
        /// <param name="flag">
        /// The flag.
        /// </param>
        public void flagEdge(int flag)
        {
            this.flagged[flag] = true;
            if (this.sibling != null)
            {
                this.sibling.flagged[flag] = true;
            }
        }

        /**
         * Sets this halfedge's <tt>flag</tt> to false. Does not change the
         * sibling's flag.
         * 
         * @param flag
         */

        /// <summary>
        /// The unflag.
        /// </summary>
        /// <param name="flag">
        /// The flag.
        /// </param>
        public void unflag(int flag)
        {
            this.flagged[flag] = false;
        }
    }

    /**
     * Represents a 2D triangulation using points (vertices) and halfedges.
     * Triangulation is performed by an interactive constrained Delaunay
     * algorithm. Point insertion uses Lawson's algorithm. Constraint insertion
     * and removal, as well as vertex removal and relocation are supported.
     *
     * For full details, see the "Interactive Constrained Delaunay Triangulation"
     * section of:
     *
     * Howison, M. CAD Tools for Creating Space-filling 3D Escher Tiles.
     * MasterвЂ™s thesis, U.C. Berkeley, Berkeley, CA, May 2009.
     * (also Tech Report EECS-2009-56,
     *     http://www.eecs.berkeley.edu/Pubs/TechRpts/2009/EECS-2009-56.html)
     *
     * Or contact the author: mark.howison <at> gmail <dot> com
     *
     * @author Mark Howison
     */

    /// <summary>
    /// The mesh.
    /// </summary>
    public class Mesh
    {
        /* drawing options */

        /// <summary>
        /// The dra w_ p t_ full.
        /// </summary>
        protected static bool DRAW_PT_FULL = true;

        /// <summary>
        /// The dra w_ p t_ indices.
        /// </summary>
        protected static bool DRAW_PT_INDICES;

        /// <summary>
        /// The dra w_ h e_ indices.
        /// </summary>
        protected static bool DRAW_HE_INDICES;

        /// <summary>
        /// The dra w_ halfedges.
        /// </summary>
        protected static bool DRAW_HALFEDGES;

        /* typical size of polygons, for initializing System.Collections.Generic.Lists */

        /// <summary>
        /// The typica l_ polygo n_ size.
        /// </summary>
        private static int TYPICAL_POLYGON_SIZE = 16;

        /* errors */

        /// <summary>
        /// The e_ exhausted.
        /// </summary>
        protected static string E_EXHAUSTED = "Exhausted halfedges or points!";

        /// <summary>
        /// The e_ missing.
        /// </summary>
        protected static string E_MISSING = "Missing halfedge or point!";

        /// <summary>
        /// The e_ identical.
        /// </summary>
        protected static string E_IDENTICAL = "Identical halfedges or points!";

        /// <summary>
        /// The e_ coincident.
        /// </summary>
        protected static string E_COINCIDENT = "Coincident points!";

        /// <summary>
        /// The e_ type.
        /// </summary>
        protected static string E_TYPE = "Incorrect type!";

        /// <summary>
        /// The e_ polygon.
        /// </summary>
        protected static string E_POLYGON = "Illegal polygonal region!";

        /// <summary>
        /// The e_ halfedge.
        /// </summary>
        protected static string E_HALFEDGE = "Mismatched halfedge!";

        /* stdout/debug/test flags */

        /// <summary>
        /// The messages.
        /// </summary>
        public static bool MESSAGES;

        /// <summary>
        /// The debug.
        /// </summary>
        public static bool DEBUG = true;

        /// <summary>
        /// The test.
        /// </summary>
        public static bool TEST;

        /* debug */

        /// <summary>
        /// The incircle error bound.
        /// </summary>
        private readonly double incircleErrorBound;

        /// <summary>
        /// The orient error bound.
        /// </summary>
        private readonly double orientErrorBound;

        /*
        // mesh data
        
        protected LinkedList<Point> points = new LinkedList<Point>();
        protected LinkedList<HalfEdge> halfEdges = new LinkedList<HalfEdge>();
        protected int nBoundary;

        // queues
        protected LinkedList<HalfEdge> delaunayQueue = new LinkedList<HalfEdge>();
        protected LinkedList<Point> removedConstraints = new LinkedList<Point>();
        protected LinkedList<Point> deleteQueue = new LinkedList<Point>();
        */

        // mesh data

        // queues
        /// <summary>
        /// The delaunay queue.
        /// </summary>
        protected List<HalfEdge> delaunayQueue = new List<HalfEdge>();

        /// <summary>
        /// The delete queue.
        /// </summary>
        protected List<Point> deleteQueue = new List<Point>();

        /// <summary>
        /// The epsilon.
        /// </summary>
        public double epsilon;

        /// <summary>
        /// The half edges.
        /// </summary>
        protected List<HalfEdge> halfEdges = new List<HalfEdge>();

        /// <summary>
        /// The n boundary.
        /// </summary>
        protected int nBoundary;

        /// <summary>
        /// The points.
        /// </summary>
        protected List<Point> points = new List<Point>();

        /// <summary>
        /// The remove constraint peg.
        /// </summary>
        private Point removeConstraintPeg;

        /// <summary>
        /// The removed constraints.
        /// </summary>
        protected List<Point> removedConstraints = new List<Point>();

        /// <summary>
        /// The toggle debug interactive.
        /// </summary>
        public bool toggleDebugInteractive;

        /// <summary>
        /// The toggle debug suspend.
        /// </summary>
        public bool toggleDebugSuspend;

        /// <summary>
        /// The toggle debug visual.
        /// </summary>
        public bool toggleDebugVisual;

        /* name and colors */
        /*
        protected String name = new String("(unnamed mesh)");
        protected Color highlightColor = Color.Red;
        protected Color warnColor = new Color(255, 240, 240);
        protected Color grayColor = new Color(127, 127, 127, 127);
        */
        /******************************************************************************/
        /* Constructors etc.
        /******************************************************************************/

        /// <summary>
        /// Initializes a new instance of the <see cref="Mesh"/> class.
        /// </summary>
        /// <param name="epsilon">
        /// The epsilon.
        /// </param>
        public Mesh(double epsilon)
        {
            this.epsilon = epsilon;
            double e = this.initMachineEpsilon();
            this.orientErrorBound = (3.0 + 16.0 * e) * e;
            this.incircleErrorBound = (10.0 + 96.0 * e) * e;

            // if (DEBUG) debugObjects = new HashSet();
        }

        /*
         * Find the machine epsilon, which is used in initializing the floating
         * point filter, the threshold at which predicate calculations become
         * unreliable and exact arithmetic is neccessary.
         *
         * This routine is adapted from Jonathan Shewchuk's predicates.c
         * exactinit() routine.
         *
         * @return  machine epsilon
         */

        /// <summary>
        /// The init machine epsilon.
        /// </summary>
        /// <returns>
        /// The init machine epsilon.
        /// </returns>
        private double initMachineEpsilon()
        {
            double half, e, check, lastcheck;

            half = 0.5;
            e = 1.0;
            check = 1.0;

            /* Repeatedly divide `epsilon' by two until it is too small to Add to    */
            /*   one without causing roundoff.  (Also check if the sum is equal to   */
            /*   the previous sum, for machines that round up instead of using exact */
            /*   rounding.  Not that this library will work on such machines anyway. */
            do
            {
                lastcheck = check;
                e *= half;
                check = 1.0 + e;
            }
            while ((check != 1.0) && (check != lastcheck));
            return e;
        }

        /******************************************************************************/
        /* Accessors & Mutators
        /******************************************************************************/

        /// <summary>
        /// The size.
        /// </summary>
        /// <returns>
        /// The size.
        /// </returns>
        public int size()
        {
            return this.points.Count;
        }

        /// <summary>
        /// The boundary size.
        /// </summary>
        /// <returns>
        /// The boundary size.
        /// </returns>
        public int boundarySize()
        {
            return this.nBoundary;
        }

        /// <summary>
        /// The clear.
        /// </summary>
        public void clear()
        {
            this.points.Clear();
            this.halfEdges.Clear();
            this.delaunayQueue.Clear();
        }

        /// <summary>
        /// The clear flags.
        /// </summary>
        /// <param name="flag">
        /// The flag.
        /// </param>
        public void clearFlags(int flag)
        {
            foreach (HalfEdge he in this.halfEdges)
            {
                he.unflag(flag);
            }
        }

        /// <summary>
        /// The contains.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <returns>
        /// The contains.
        /// </returns>
        public bool contains(Point p)
        {
            return this.points.Contains(p);
        }

        /// <summary>
        /// The index of.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <returns>
        /// The index of.
        /// </returns>
        public int IndexOf(Point p)
        {
            return this.points.IndexOf(p);
        }

        /// <summary>
        /// The get point.
        /// </summary>
        /// <param name="i">
        /// The i.
        /// </param>
        /// <returns>
        /// </returns>
        public Point getPoint(int i)
        {
            // return points[i);
            return this.points[i];
        }

        /*
        public Point3d[] getPointCoordinates3d()
        {
            Point3d[] coords = new Point3d[points.Count];
            int i=0;
            foreach (Point p in points) {
                coords[i++] = p.coords3d;
            }
            return coords;
        }
        */

        /// <summary>
        /// The get face coordinates.
        /// </summary>
        /// <returns>
        /// </returns>
        public Point2d[] getFaceCoordinates()
        {
            var facePoints = new List<Point2d>();

            // reset 'used' flags
            this.clearFlags(HalfEdge.FLAG_READ);

            // find the faces
            foreach (HalfEdge he0 in this.halfEdges)
            {
                if (he0.isFlagged(HalfEdge.FLAG_READ))
                {
                    continue;
                }

                HalfEdge he1 = he0.next;
                HalfEdge he2 = he1.next;

                // the face is oriented CCW
                facePoints.Add(new Point2d(he0.origin));
                facePoints.Add(new Point2d(he1.origin));
                facePoints.Add(new Point2d(he2.origin));

                // mark these half edges as used
                he0.flag(HalfEdge.FLAG_READ);
                he1.flag(HalfEdge.FLAG_READ);
                he2.flag(HalfEdge.FLAG_READ);
            }

            return facePoints.ToArray();
        }

        /*
            public Point3d[] getFaceCoordinates3d()
            {
                System.Collections.Generic.List<Point3d> coords = new System.Collections.Generic.List();
                // reset 'used' flags
                clearFlags(HalfEdge.FLAG_READ);
                // find the faces
                foreach (HalfEdge he0 in halfEdges) {
                    if (he0.isFlagged(HalfEdge.FLAG_READ)) continue;
                    HalfEdge he1 = he0.next;
                    HalfEdge he2 = he1.next;
                    // the face is oriented CCW
                    coords.Add(new Point3d(he0.origin.coords3d));
                    coords.Add(new Point3d(he1.origin.coords3d));
                    coords.Add(new Point3d(he2.origin.coords3d));
                    // mark these half edges as used
                    he0.flag(HalfEdge.FLAG_READ);
                    he1.flag(HalfEdge.FLAG_READ);
                    he2.flag(HalfEdge.FLAG_READ);
                }
                //// assert coords.Count % 3 == 0 : error("Illegal face list!");
                return coords.toArray(new Point3d[0]);
            }
            */

        /// <summary>
        /// The get edges.
        /// </summary>
        /// <returns>
        /// </returns>
        public Edge[] getEdges()
        {
            var edges = new List<Edge>();
            this.clearFlags(HalfEdge.FLAG_READ);

            // find the faces
            foreach (HalfEdge he in this.halfEdges)
            {
                if (he.isFlagged(HalfEdge.FLAG_READ))
                {
                    continue;
                }

                var e = new Edge(new Point(he.origin), new Point(he.next.origin));
                e.type = he.getType();
                edges.Add(e);
                he.flag(HalfEdge.FLAG_READ);
            }

            return edges.ToArray();
        }

        /*
                public void setName(String s)
                {
                    name = s;
                }
                */
        /*
        public void setDebugFrame(JFrame f)
        {
            debugFrame = f;
        }

        public void setDebugPanel(JPanel p)
        {
            debugPanel = p;
        }
        */

        /// <summary>
        /// The translate.
        /// </summary>
        /// <param name="x">
        /// The x.
        /// </param>
        /// <param name="y">
        /// The y.
        /// </param>
        public void translate(double x, double y)
        {
            foreach (Point p in this.points)
            {
                p.X += x;
                p.Y += y;
            }
        }

        /******************************************************************************/
        /* Insertion methods
        /******************************************************************************/

        /// <summary>
        /// The link boundary.
        /// </summary>
        /// <param name="pts">
        /// The pts.
        /// </param>
        public static void linkBoundary(BPoint[] pts)
        {
            int s = pts.Length;
            for (int i = 0; i < s; i++)
            {
                pts[i].next = pts[(i + 1) % s];
                pts[i].prev = pts[(i - 1 + s) % s];
            }
        }

        // expects a list of points in *clockwise* order for specifiying
        // an initial boundary for the mesh
        /// <summary>
        /// The init.
        /// </summary>
        /// <param name="pts">
        /// The pts.
        /// </param>
        public void init(Point[] pts)
        {
            int s = pts.Length;

            //// assert s >= 3 : error("Initialization requires at least 3 points!");
            var polygon = new List<HalfEdge>(s);
            this.clear();
            foreach (Point p in pts)
            {
                this.points.Add(p);
                p.setType(Point.BOUNDARY);
                p.he = new HalfEdge(p, HalfEdge.BOUNDARY);
                this.halfEdges.Add(p.he);
            }

            for (int i = 0; i < s; i++)
            {
                this.halfEdges[i].next = this.halfEdges[(s + i - 1) % s];
            }

            this.nBoundary = s;
            for (int i = 0; i < s; i++)
            {
                polygon.Add(this.halfEdges[s - 1 - i]);
            }

            fillGeneralPolygon(polygon);
            this.delaunayQueue.AddRange(polygon);

            /*
            foreach(HalfEdge edge in polygon){
                delaunayQueue.AddLast(edge);
            }
             */
            this.updateDelaunay();
        }

        /// <summary>
        /// The add half edge.
        /// </summary>
        /// <param name="origin">
        /// The origin.
        /// </param>
        /// <param name="destination">
        /// The destination.
        /// </param>
        /// <returns>
        /// </returns>
        private HalfEdge addHalfEdge(Point origin, Point destination)
        {
            var he1 = new HalfEdge(origin);
            var he2 = new HalfEdge(destination);
            he1.sibling = he2;
            he2.sibling = he1;
            this.halfEdges.Add(he1);
            this.halfEdges.Add(he2);
            if (MESSAGES)
            {
                this.message("Added halfedges %d and %d.", this.halfEdges.IndexOf(he1), this.halfEdges.IndexOf(he2));
            }

            return he1;
        }

        /*
         * Adds an edge connecting the origins of he1 and he2.
         *
         * @return  the halfedge from he1.origin to he2.origin
         */

        /// <summary>
        /// The add edge.
        /// </summary>
        /// <param name="he1">
        /// The he 1.
        /// </param>
        /// <param name="he2">
        /// The he 2.
        /// </param>
        /// <param name="he1prev">
        /// The he 1 prev.
        /// </param>
        /// <param name="he2prev">
        /// The he 2 prev.
        /// </param>
        /// <returns>
        /// </returns>
        private HalfEdge addEdge(HalfEdge he1, HalfEdge he2, HalfEdge he1prev, HalfEdge he2prev)
        {
            //// assert he1prev.next == he1 : error(E_HALFEDGE,he1,he1prev);
            //// assert he2prev.next == he2 : error(E_HALFEDGE,he2,he2prev);
            //// assert he1 != he2 : error(E_COINCIDENT);
            //// assert he1.origin != he2.origin: error(E_COINCIDENT);
            HalfEdge heAdd = this.addHalfEdge(he1.origin, he2.origin);
            this.delaunayQueue.Add(heAdd);
            heAdd.next = he2;
            he1prev.next = heAdd;
            heAdd.sibling.next = he1;
            he2prev.next = heAdd.sibling;
            return heAdd;
        }

        /// <summary>
        /// The add boundary point.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <param name="he0">
        /// The he 0.
        /// </param>
        /// <returns>
        /// </returns>
        public Point addBoundaryPoint(Point p, HalfEdge he0)
        {
            //// assert halfEdges.contains(he0) : error(E_MISSING);
            //// assert between(he0.origin,he0.next.origin,p) :
            this.error("Adding boundary point that doesn't lie on boundary!");

            HalfEdge he1, he2, he3;

            /* check for coincidence with the endpoints of he0 */
            if (this.coincident(p, he0.origin))
            {
                /*
                if (MESSAGES) message(
                        "Boundary point is within epsilon of %d.",
                        points.IndexOf(he0.origin));
                 */
                return he0.origin;
            }

            if (this.coincident(p, he0.next.origin))
            {
                /*
                if (MESSAGES) message(
                        "Boundary point is within epsilon of %d.",
                        points.IndexOf(he0.next.origin));*/
                return he0.next.origin;
            }

            this.points.Add(p);

            // if (MESSAGES) message("Adding boundary point %d.", points.IndexOf(p));
            he2 = he0.next;
            he3 = he2.next;

            /* split the existing boundary */
            he1 = new HalfEdge(p, HalfEdge.BOUNDARY);
            this.halfEdges.Add(he1);

            /* link halfedges */
            p.he = he1;
            he0.next = he1;
            he1.next = he2;
            this.fillQuadrilateral(he0);
            this.updateDelaunay();
            this.nBoundary++;

            // if (TEST) test();
            return p;
        }

        /// <summary>
        /// Нахождение квадрата расстояния между двумя точками
        /// </summary>
        /// <param name="p1">
        /// Точка 1
        /// </param>
        /// <param name="p2">
        /// Точка 2
        /// </param>
        /// <returns>
        /// расстояние между точками
        /// </returns>
        public static double FindSquareDist(Point p1, Point p2)
        {
            double Dist = ((p1.X - p2.X) * (p1.X - p2.X)) + ((p1.Y - p2.Y) * (p1.Y - p2.Y));
            return Dist;
        }

        /// <summary>
        /// The add interior point.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <returns>
        /// </returns>
        public Point addInteriorPoint(Point p)
        {
            double dist, min;
            Point pNearest = null;
            FaceWalk walk;

            // if (DEBUG) debugView(p, "addInteriorPoint: p");
            /* find the closest point to p */
            min = Double.MaxValue;
            foreach (Point pTest in this.points)
            {
                dist = FindSquareDist(p, pTest);

                /* abort if the point is within epsilon of an existing point */
                if (dist < this.epsilon)
                {
                    /*
                    if (MESSAGES) message(
                            "Point is within epsilon of %d.",
                            points.IndexOf(pTest));
                    if (DEBUG) debugView(pTest, "addInteriorPoint: within epsilon");*/
                    return pTest;
                }
                else if (dist < min)
                {
                    min = dist;
                    pNearest = pTest;
                }
            }

            // if (MESSAGES) message("Closest point is ", points.IndexOf(pNearest));
            // if (DEBUG) debugView(pNearest, "addInteriorPoint: closest point");
            /* find face containing p, starting at pNearest */
            walk = this.findFace(pNearest.he, p);
            if (MESSAGES)
            {
                this.message("Point is within face with halfedge %d.", this.halfEdges.IndexOf(walk.he));
            }

            this.points.Add(p);
            if (walk.status == FaceWalk.COINCIDENT)
            {
                this.splitEdge(p, walk.he);
            }
            else
            {
                this.splitFace(p, walk.he);
            }

            this.updateDelaunay();
            return p;
        }

        /// <summary>
        /// The add constraint.
        /// </summary>
        /// <param name="pStart">
        /// The p start.
        /// </param>
        /// <param name="pEnd">
        /// The p end.
        /// </param>
        /// <returns>
        /// The add constraint.
        /// </returns>
        public bool addConstraint(Point pStart, Point pEnd)
        {
            //// assert points.contains(pStart) : error(E_MISSING);
            //// assert points.contains(pEnd) : error(E_MISSING);
            //// assert pStart != pEnd : error("Identical points!");
            //// assert ! coincident(pStart,pEnd) : error("Coincident points!");
            int i;
            Point pSearch0, pSearch1;
            HalfEdge heSearch, heStart, heStartPrev;
            FaceWalk walk;

            if (MESSAGES)
            {
                this.message("Constraining %d -> %d.", this.points.IndexOf(pStart), this.points.IndexOf(pEnd));
            }

            /* find the halfedge at pStart that lies on or below the constraint */
            walk = this.startFaceWalk(pStart, pEnd);

            /* check for trivial condition where the edge already exists */
            if (walk.status == FaceWalk.COINCIDENT)
            {
                return this.constrainEdge(walk.he);
            }

            /* clear edges that intersect the constraint */
            heStart = walk.he;
            heStartPrev = this.findPrevious(heStart);
            heSearch = heStart.next;
            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                pSearch0 = heSearch.origin;
                pSearch1 = heSearch.next.origin;

                // if (DEBUG) debugView(heSearch, "addConstraint: walking");
                /* check for termination */
                if (pSearch1 == pEnd)
                {
                    if (MESSAGES)
                    {
                        this.message("Found constraint end at halfedge %d.", this.halfEdges.IndexOf(heSearch));
                    }

                    break;
                }

                //// assert ! coincident(pSearch1,pStart) :
                // error(E_COINCIDENT, pSearch1, pStart);
                //// assert ! coincident(pSearch1,pEnd) :
                // error(E_COINCIDENT, pSearch1, pEnd);
                /* check for collinearity */
                // if (between(pStart,pEnd,pSearch1)) {
                /* split the constraint in two with pSearch1 as the midpoint */
                /*addConstraintEdge(heStart,heSearch.next,heStartPrev,heSearch);
                return addConstraint(pSearch1,pEnd);
            }*/
                /* check for intersection */
                if (this.intersect(pStart, pEnd, pSearch0, pSearch1))
                {
                    //// assert ! heSearch.isType(HalfEdge.BOUNDARY) :
                    this.error("Constraint crosses boundary edge!");
                    if (heSearch.isType(HalfEdge.AUXILARY))
                    {
                        this.removeEdge(heSearch);
                        heSearch = heSearch.sibling;
                    }
                    else if (heSearch.isType(HalfEdge.CONSTRAINT))
                    {
                        if (MESSAGES)
                        {
                            this.message("Constraint-constraint intersection found.");
                        }

                        Point2d p = intersection(pStart, pEnd, pSearch0, pSearch1);
                        this.splitConstraint(heSearch, p);
                        this.addConstraintEdge(heStart, heSearch.next, heStartPrev, heSearch);

                        /* reset the starting point */
                        heSearch = heSearch.sibling;
                        heStart = heSearch;
                        heStartPrev = this.findPrevious(heStart);
                        pStart = heSearch.origin;
                    }
                }

                heSearch = heSearch.next;
            }

            //// assert i < halfEdges.Count : error(E_EXHAUSTED,pStart,pEnd);
            this.addConstraintEdge(heStart, heSearch.next, heStartPrev, heSearch);
            this.updateDelaunay();
            return true;
        }

        /// <summary>
        /// The add constraint edge.
        /// </summary>
        /// <param name="he1">
        /// The he 1.
        /// </param>
        /// <param name="he2">
        /// The he 2.
        /// </param>
        /// <param name="he1prev">
        /// The he 1 prev.
        /// </param>
        /// <param name="he2prev">
        /// The he 2 prev.
        /// </param>
        private void addConstraintEdge(HalfEdge he1, HalfEdge he2, HalfEdge he1prev, HalfEdge he2prev)
        {
            //// assert(he1prev.next == he1, error(E_HALFEDGE));
            //// assert he2prev.next == he2 : error(E_HALFEDGE);
            HalfEdge heAdd = this.addEdge(he1, he2, he1prev, he2prev);
            this.constrainEdge(heAdd);
            this.fillEdgeVisiblePolygon(heAdd);
            this.fillEdgeVisiblePolygon(heAdd.sibling);
        }

        /******************************************************************************/
        /* Polygon filling methods
        /******************************************************************************/

        /// <summary>
        /// The fill quadrilateral.
        /// </summary>
        /// <param name="he1">
        /// The he 1.
        /// </param>
        protected void fillQuadrilateral(HalfEdge he1)
        {
            HalfEdge he2, he3, he4;

            he2 = he1.next;
            he3 = he2.next;
            he4 = he3.next;

            //// assert he4.next == he1;
            if (this.orient(he1.origin, he3.origin, he2.origin) * this.orient(he1.origin, he3.origin, he4.origin) < 0)
            {
                this.addEdge(he1, he3, he4, he2);
            }
            else
            {
                this.addEdge(he2, he4, he1, he3);
            }
        }

        /// <summary>
        /// The construct polygon.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        /// <returns>
        /// </returns>
        protected List<HalfEdge> constructPolygon(HalfEdge he)
        {
            //// assert halfEdges.contains(he) : error(E_MISSING);
            int i;
            HalfEdge heSearch;
            var polygon = new List<HalfEdge>(TYPICAL_POLYGON_SIZE);

            polygon.Add(he);
            heSearch = he.next;
            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                if (heSearch == he)
                {
                    break;
                }

                polygon.Add(heSearch);

                // if (DEBUG) debugView(heSearch, "constructPolygon: next");
                heSearch = heSearch.next;
            }

            //// assert i < halfEdges.Count : error(E_EXHAUSTED,he,he.origin);
            return polygon;
        }

        /// <summary>
        /// The fill general polygon.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        protected void fillGeneralPolygon(HalfEdge he)
        {
            this.fillGeneralPolygon(this.constructPolygon(he));
        }

        /// <summary>
        /// The fill general polygon.
        /// </summary>
        /// <param name="polygon">
        /// The polygon.
        /// </param>
        protected void fillGeneralPolygon(List<HalfEdge> polygon)
        {
            this.fillGeneralPolygonRecurse(polygon);
            this.delaunayQueue.AddRange(polygon);
        }

        /// <summary>
        /// The fill general polygon recurse.
        /// </summary>
        /// <param name="polygon">
        /// The polygon.
        /// </param>
        private void fillGeneralPolygonRecurse(List<HalfEdge> polygon)
        {
            //// assert polygon.Count >= 3 : error("Illegal size!");
            int n, s;
            Point p0, p1, p2;
            HalfEdge heTest0, heTest1, heAdd;

            /* Assumes a Jordan (simple) polygon with n>3 sides!
             * (i.e. no sides intersect) */
            if (MESSAGES)
            {
                this.message("Filling polygon with %d sides.");
            }

            // if (DEBUG) debugView(polygon[0), "fillPolygon: start");
            s = polygon.Count;
            if (s > 3)
            {
                /* A Jordan polygon always has two non-overlapping ears.
                 * We iterate over all possible ear edges,
                 * i.e. those between vertices i and i+2 in the polygon. */
                p0 = p1 = p2 = null;
                heTest0 = null;
                n = 0;

                for (int i = 0; i < s; i++)
                {
                    bool wasBreak = false;
                    n = i;
                    heTest0 = polygon[i];
                    p0 = heTest0.origin;
                    p1 = heTest0.next.origin;
                    p2 = polygon[(i + 2) % s].origin;

                    // check that the ear edge p0->p2 lies strictly
                    // inside the polygon, i.e. to the left of p0->p1
                    if (this.orient(p0, p1, p2) > 0 && (!this.between(p0, p2, p1)))
                    {
                        // check for intersections or points that lie too
                        // close to the ear edge
                        heTest1 = heTest0.next.next;
                        for (int j = 0; j < (s - 3); j++)
                        {
                            if (this.intersectProper(p0, p2, heTest1.origin, heTest1.next.origin))
                            {
                                // continue edgeWalk;
                                wasBreak = true;
                                break;
                            }

                            heTest1 = heTest1.next;
                        }

                        if (wasBreak)
                        {
                            continue;
                        }

                        break;
                    }
                }

                heAdd = this.addHalfEdge(p0, p2);
                this.delaunayQueue.Add(heAdd);

                // link halfedges in the ear
                heAdd.sibling.next = heTest0;
                polygon[(n + 1) % s].next = heAdd.sibling;

                // link halfedges in the remaining polygon of size s-1
                heAdd.next = polygon[(n + 2) % s];
                polygon[(n + s - 1) % s].next = heAdd;
                if (s > 4)
                {
                    var polygon0 = new List<HalfEdge>();
                    for (int j = 0; j < (s - 1); j++)
                    {
                        polygon0.Add(heAdd);
                        heAdd = heAdd.next;
                    }

                    this.fillGeneralPolygonRecurse(polygon0);
                }
            }

            // if (TEST) test();
        }

        /// <summary>
        /// The fill edge visible polygon.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        protected void fillEdgeVisiblePolygon(HalfEdge he)
        {
            List<HalfEdge> polygon = this.constructPolygon(he);
            this.fillEdgeVisiblePolygonRecurse(polygon);
            this.delaunayQueue.AddRange(polygon);
        }

        /// <summary>
        /// The fill edge visible polygon recurse.
        /// </summary>
        /// <param name="polygon">
        /// The polygon.
        /// </param>
        private void fillEdgeVisiblePolygonRecurse(List<HalfEdge> polygon)
        {
            //// assert polygon.Count >= 3 : error("Illegal size!");
            int i, c, s;
            Point pa, pb, pc;
            HalfEdge heAdd;

            if (MESSAGES)
            {
                this.message("Filling polygon with %d sides.", polygon.Count);
            }

            // if (DEBUG) debugView(polygon[0], "fillPolygon: start");
            s = polygon.Count;
            if (s > 3)
            {
                pa = polygon[0].origin;
                pb = polygon[1].origin;
                pc = polygon[2].origin;
                c = 2;
                for (i = 3; i < s; i++)
                {
                    Point p = polygon[i].origin;
                    if (this.incircle(pa, pb, pc, p) > 0)
                    {
                        pc = p;
                        c = i;
                    }
                }

                // if (DEBUG) debugView(polygon[c],                      "fillPolygon: largest circumcircle");
                /* Add edge pa -> pc */
                if (c < (s - 1))
                {
                    heAdd = this.addEdge(polygon[0], polygon[c], polygon[s - 1], polygon[c - 1]);
                    this.fillEdgeVisiblePolygonRecurse(this.constructPolygon(heAdd));
                }

                /* Add edge pb -> pc */
                if (c > 2)
                {
                    heAdd = this.addEdge(polygon[1], polygon[c - 1].next, polygon[0], polygon[c - 1]);
                    this.fillEdgeVisiblePolygonRecurse(this.constructPolygon(heAdd.sibling));
                }
            }

            // if (TEST) test();
        }

        /******************************************************************************/
        /* Splitting methods
        /******************************************************************************/

        /*
         * Splits the two faces sharing edge he into four faces by inserting
         * point p.
         */

        /// <summary>
        /// The split edge.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <param name="he">
        /// The he.
        /// </param>
        private void splitEdge(Point p, HalfEdge he)
        {
            HalfEdge he1, he2, he3;
            HalfEdge heAdd1, heAdd2, heAdd3;

            if (MESSAGES)
            {
                this.message("Splitting edge %d.", this.halfEdges.IndexOf(he));
            }

            //// assert ! he.isType(HalfEdge.BOUNDARY) :
            this.error("Attempting to split boundary edge!");

            // he1 = he.next.next;
            he1 = this.findPrevious(he);
            he2 = he.sibling.next;

            // he3 = he2.next;
            he3 = this.findPrevious(he.sibling);

            // split the halfedge
            he.origin = p;
            p.he = he;

            // Add halfedges
            heAdd1 = this.addHalfEdge(p, he1.origin);
            heAdd2 = this.addHalfEdge(p, he2.origin);
            heAdd3 = this.addHalfEdge(p, he3.origin);

            // link halfedges
            heAdd1.next = he1;

            // heAdd2.next = he2;
            // heAdd3.next = he3;
            he.next.next = heAdd1.sibling;
            he1.next = heAdd2.sibling;
            he2.next = heAdd3.sibling;
            heAdd1.sibling.next = he;
            heAdd2.sibling.next = heAdd1;
            heAdd3.sibling.next = heAdd2;
            he.sibling.next = heAdd3;

            // update the point->halfedge pointers
            this.updateHalfEdge(he2);

            // Add halfedges to delaunay test
            this.delaunayQueue.Add(he);
            this.delaunayQueue.Add(he1);
            this.delaunayQueue.Add(he2);
            this.delaunayQueue.Add(he3);
        }

        /*
         * Insert point p into the face with halfedge he1, splitting it into
         * three faces.
         */

        /// <summary>
        /// The split face.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <param name="he1">
        /// The he 1.
        /// </param>
        private void splitFace(Point p, HalfEdge he1)
        {
            //// assert halfEdges.contains(he1) : error(E_MISSING);
            HalfEdge he2, he3;
            HalfEdge heAdd1, heAdd2, heAdd3;

            if (MESSAGES)
            {
                this.message("Adding interior point inside face.");
            }

            he2 = he1.next;
            he3 = he2.next;

            /* Add new halfedges */
            heAdd1 = this.addHalfEdge(p, he1.origin);
            heAdd2 = this.addHalfEdge(p, he2.origin);
            heAdd3 = this.addHalfEdge(p, he3.origin);

            /* link half edges */
            p.he = heAdd1;
            heAdd1.next = he1;
            heAdd2.next = he2;
            heAdd3.next = he3;
            he1.next = heAdd2.sibling;
            he2.next = heAdd3.sibling;
            he3.next = heAdd1.sibling;
            heAdd1.sibling.next = heAdd3;
            heAdd3.sibling.next = heAdd2;
            heAdd2.sibling.next = heAdd1;

            /* Add halfedges to delaunay test */
            this.delaunayQueue.Add(heAdd1);
            this.delaunayQueue.Add(heAdd2);
            this.delaunayQueue.Add(heAdd3);
            this.delaunayQueue.Add(he1);
            this.delaunayQueue.Add(he2);
            this.delaunayQueue.Add(he3);
        }

        /// <summary>
        /// The split constraint.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        /// <param name="p">
        /// The p.
        /// </param>
        private void splitConstraint(HalfEdge he, Point2d p)
        {
            int i;
            Point p0;
            HalfEdge he0, heTest;

            if (MESSAGES)
            {
                this.message("Splitting constraint edge %d.", this.halfEdges.IndexOf(he));
            }

            //// assert ! he.isType(HalfEdge.BOUNDARY) :
            this.error("Attempting to split a boundary edge!");

            // Add point
            p0 = new Point(p);
            this.points.Add(p0);
            if (MESSAGES)
            {
                this.message("Adding constraint intersection point %d.", this.points.IndexOf(p0));
            }

            // Add halfedge
            he0 = this.addHalfEdge(p0, he.sibling.origin);
            he0.constrain();
            he0.sibling.constrain();
            he.sibling.origin = p0;

            // update point->halfedge pointers
            p0.he = he0;
            this.updateHalfEdge(he0.sibling);

            // link halfedges
            he0.next = he.next;
            he.next = he0;
            he0.sibling.next = he.sibling;

            // find halfedge pointing to he.sibling
            heTest = he.sibling;
            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                if (heTest.next == he.sibling)
                {
                    heTest.next = he0.sibling;
                    break;
                }

                heTest = heTest.next;
            }

            //// assert i < halfEdges.Count : error(E_EXHAUSTED);
        }

        /******************************************************************************/
        /* Removal methods
        /******************************************************************************/

        /// <summary>
        /// The remove boundary point.
        /// </summary>
        /// <param name="bp">
        /// The bp.
        /// </param>
        public void removeBoundaryPoint(BPoint bp)
        {
            BPoint bp1, bp2;

            /* The boundary is oriented clockwise, opposite of the halfedge
             * orientation (counter-clockwise). Therefore, the "previous" point
             * to bp is actually the next point along the boundary. */
            this.removeBoundaryPoint(bp, bp.next);

            /* relink boundary */
            bp1 = bp.prev;
            bp2 = bp.next;
            bp1.next = bp2;
            bp2.prev = bp1;
        }

        /// <summary>
        /// The remove boundary point.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <param name="pPrev">
        /// The p prev.
        /// </param>
        private void removeBoundaryPoint(Point p, Point pPrev)
        {
            int i;
            Point pNext;
            HalfEdge he;

            if (MESSAGES)
            {
                this.message("Removing boundary point %d.", this.points.IndexOf(p));
            }

            pNext = p.he.next.origin;
            he = p.he.getPrev();
            this.deleteQueue.Clear();
            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                if (he.isType(HalfEdge.BOUNDARY))
                {
                    break;
                }

                if (this.orient(pPrev, pNext, he.origin) <= 0)
                {
                    this.deleteQueue.Add(he.origin);
                }

                this.removeEdge(he);
                he = p.he.getPrev();
            }

            //// assert i < halfEdges.Count : error(E_EXHAUSTED);
            if (MESSAGES)
            {
                this.message("Performing flood delete to Remove interior points.");
            }

            var bounds = new Point[3];
            bounds[0] = pPrev;
            bounds[1] = pNext;
            bounds[2] = p;
            this.floodDelete(bounds);
            pPrev.he.next = p.he.next;
            this.halfEdges.Remove(p.he);
            p.he = null;
            this.points.Remove(p);
            p.setType(Point.DELETED);

            /* stitch the polygon back together */
            this.fillEdgeVisiblePolygon(pPrev.he);
            this.updateDelaunay();
            this.nBoundary--;
        }

        /// <summary>
        /// The remove interior point.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        public void removeInteriorPoint(Point p)
        {
            //// assert points.contains(p) : error(E_MISSING);
            int i;
            Point p1, p2, p3;
            HalfEdge heSearch, heFlip;
            var star = new List<HalfEdge>();

            if (MESSAGES)
            {
                this.message("Removing interior point.");
            }

            /* construct the star of halfedges around p */
            heSearch = p.he;
            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                star.Add(heSearch);
                heSearch = heSearch.next.next.sibling;
                if (heSearch == p.he)
                {
                    break;
                }
            }

            //// assert i < halfEdges.Count : error(E_EXHAUSTED);

            //// assert star.Count >= 3;
            if (star.Count == 3)
            {
                foreach (HalfEdge he in star)
                {
                    this.removeEdge(he);
                }
            }
            else
            {
                while (star.Count > 4)
                {
                    heFlip = star[0];
                    star.RemoveAt(0);
                    p1 = heFlip.sibling.origin;
                    p2 = heFlip.next.next.origin;
                    p3 = heFlip.sibling.next.next.origin;
                    if (this.orient(p2, p3, p) * this.orient(p2, p3, p1) < 0)
                    {
                        this.flipEdge(heFlip);
                    }
                    else
                    {
                        star.Add(heFlip);
                    }
                }

                //// assert star.Count == 4;
                heSearch = star[0].next;
                foreach (HalfEdge he in star)
                {
                    this.removeEdge(he);
                }

                this.fillQuadrilateral(heSearch);
            }

            p.he = null;
            this.points.Remove(p);
            p.setType(Point.DELETED);
            this.updateDelaunay();
        }

        /*
         * Removes an interior point without refilling.
         *
         * @params p
         */

        /// <summary>
        /// The remove point.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        private void removePoint(Point p)
        {
            //// assert points.contains(p) : error(E_MISSING);
            //// assert ! p.isType(Point.DELETED) : error("Re-removing point!");
            int i;
            HalfEdge he = p.he;

            if (MESSAGES)
            {
                this.message("Removing point %d.", this.points.IndexOf(p));
            }

            // if (DEBUG) debugView(he, "removePoint: he");
            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                //// assert he.origin == p : error("Mismatched halfedge!");
                this.removeEdge(he);
                if (p.he == null)
                {
                    break;
                }

                he = he.sibling.next;
            }

            //// assert i < halfEdges.Count : error(E_EXHAUSTED);
            this.points.Remove(p);
            p.setType(Point.DELETED);
            p.he = null;
        }

        /// <summary>
        /// The remove edge.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        private void removeEdge(HalfEdge he)
        {
            //// assert halfEdges.contains(he) : error(E_MISSING);
            //// assert ! he.isType(HalfEdge.BOUNDARY) : error(E_TYPE);
            HalfEdge hePrev, heSibPrev;

            // if (DEBUG) debugView(he, "removeEdge: he");
            if (MESSAGES)
            {
                this.message("Removing edge %d.", this.halfEdges.IndexOf(he));
            }

            hePrev = this.findPrevious(he);
            heSibPrev = this.findPrevious(he.sibling);

            /* Remove halfedges */
            this.halfEdges.Remove(he);
            this.halfEdges.Remove(he.sibling);
            this.delaunayQueue.Remove(he);
            this.delaunayQueue.Remove(he.sibling);

            /* cache the constraints */
            if (he.isType(HalfEdge.CONSTRAINT))
            {
                this.removedConstraints.Add(he.next.origin);
            }

            /* update point->halfedge pointers */
            if (he.sibling == hePrev)
            {
                /* this was the last halfedge eminating from he.origin */
                // if (DEBUG) debugView(he.origin, "removeEdge: orphan point");
                he.origin.he = null;
                this.updateHalfEdge(he.next);
            }
            else if (he.next == he.sibling)
            {
                /* this was the last halfedge eminating from he.sibling.origin */
                // if (DEBUG) debugView(he.next.origin, "removeEdge: orphan point");
                he.next.origin.he = null;
                this.updateHalfEdge(he.sibling.next);
            }
            else
            {
                this.updateHalfEdge(he.next);
                this.updateHalfEdge(he.sibling.next);
            }

            /* relink halfedges */
            hePrev.next = he.sibling.next;
            heSibPrev.next = he.next;
        }

        /// <summary>
        /// The flood delete.
        /// </summary>
        /// <param name="bounds">
        /// The bounds.
        /// </param>
        private void floodDelete(Point[] bounds)
        {
            //// assert bounds.Length >= 3 : error("Illegal bounds!");
            int i;
            int[] types;
            bool inside;
            Point p1, p2;
            HalfEdge heTest;
            var boundsPath = new GraphicsPath();
            boundsPath.StartFigure();

            // GeneralPath boundsPath = new GeneralPath();
            types = new int[bounds.Length];
            for (i = 0; i < bounds.Length; i++)
            {
                types[i] = bounds[i].type;
            }

            foreach (Point p in bounds)
            {
                p.type = Point.BOUNDS;
            }

            // boundsPath.moveTo(bounds[0].X, bounds[0].Y);
            boundsPath.AddLine((float)bounds[0].X, (float)bounds[0].Y, (float)bounds[1].X, (float)bounds[1].Y);
            for (i = 2; i < bounds.Length; i++)
            {
                boundsPath.AddLine(
                    (float)bounds[i - 1].X, (float)bounds[i - 1].Y, (float)bounds[i].X, (float)bounds[i].Y);
            }

            boundsPath.CloseFigure();

            // boundsPath.closePath();
            while (this.deleteQueue.Count > 0)
            {
                p1 = this.deleteQueue[0];
                this.deleteQueue.RemoveAt(0);
                bool contain = false;
                foreach (PointF point in boundsPath.PathPoints)
                {
                    if (((p1.X - point.X) < float.Epsilon) && ((p1.Y - point.Y) < float.Epsilon))
                    {
                        contain = true;
                        break;
                    }
                }

                if (contain && !p1.isType(Point.BOUNDS))
                {
                    heTest = p1.he;
                    for (i = 0; i <= this.halfEdges.Count; i++)
                    {
                        if (p1.he == null)
                        {
                            break;
                        }

                        p2 = heTest.next.origin;
                        if (!p2.isType(Point.BOUNDS))
                        {
                            inside = true;
                            for (int j = 1; j < bounds.Length; j++)
                            {
                                if (this.intersect(p1, p2, bounds[j - 1], bounds[j]))
                                {
                                    /* check if p2 lies strictly outside the bounds
                                     * (i.e. is not colinear with the bounds) */
                                    if (this.orient(bounds[j - 1], bounds[j], p2) != 0)
                                    {
                                        inside = false;
                                    }

                                    break;
                                }
                            }

                            if (inside)
                            {
                                this.deleteQueue.Add(p2);
                            }
                        }

                        this.removeEdge(heTest);
                        heTest = heTest.sibling.next;
                    }

                    //// assert i < halfEdges.Count : error(E_EXHAUSTED);
                    this.points.Remove(p1);
                    p1.type = Point.DELETED;
                    p1.he = null;
                }
            }

            for (i = 0; i < bounds.Length; i++)
            {
                bounds[i].type = types[i];
            }
        }

        /// <summary>
        /// The clear new boundary edge.
        /// </summary>
        /// <param name="pStart">
        /// The p start.
        /// </param>
        /// <param name="pEnd">
        /// The p end.
        /// </param>
        private void clearNewBoundaryEdge(Point pStart, Point pEnd)
        {
            //// assert points.contains(pStart) : error(E_MISSING);
            //// assert points.contains(pEnd) : error(E_MISSING);
            //// assert pStart != pEnd : error(E_IDENTICAL);
            FaceWalk walk;

            if (MESSAGES)
            {
                this.message("Clearing new boundary edge.");
            }

            walk = this.startFaceWalk(pStart, pEnd);

            // check for trivial case
            if (walk.status == FaceWalk.COINCIDENT)
            {
                this.removeEdge(walk.he);
            }
            else
            {
                int i;
                Point pSearch0, pSearch1;
                HalfEdge heSearch;

                heSearch = walk.he;

                // if (DEBUG) debugView(                        heSearch,                        "clearNewBoundaryEdge: start face walk");
                for (i = 0; i <= this.halfEdges.Count; i++)
                {
                    // if (DEBUG) debugView(heSearch, "clearNewBoundaryEdge: walking");
                    pSearch0 = heSearch.origin;
                    pSearch1 = heSearch.next.origin;

                    // check for termination
                    if (pSearch1 == pEnd)
                    {
                        break;
                    }

                    // check for collinearity
                    if (this.between(pStart, pEnd, pSearch1))
                    {
                        if (pSearch1.isType(Point.BOUNDARY))
                        {
                            heSearch = heSearch.next;
                        }
                        else
                        {
                            this.deleteQueue.Add(pSearch1);
                            heSearch = heSearch.sibling.next;
                        }
                    }
                        
                        
                        // check for intersection
                    else if (this.intersectProper(pStart, pEnd, pSearch0, pSearch1))
                    {
                        this.removeEdge(heSearch);
                        this.deleteQueue.Add(pSearch0);
                        heSearch = heSearch.sibling.next;
                    }
                    else
                    {
                        heSearch = heSearch.next;
                    }
                }

                //// assert i < halfEdges.Count : error(E_EXHAUSTED);
            }
        }

        /******************************************************************************/
        /* Update methods
        /******************************************************************************/

        /// <summary>
        /// The update half edge.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        protected void updateHalfEdge(HalfEdge he)
        {
            //// assert halfEdges.contains(he) : error(E_MISSING);
            if (he.origin.isType(Point.INTERIOR))
            {
                he.origin.he = he;
            }
        }

        /// <summary>
        /// The update delaunay.
        /// </summary>
        public void updateDelaunay()
        {
            Point2d p1, p2, p3, p4;

            if (MESSAGES)
            {
                this.message("Testing Delaunay queue with %d halfedges.", this.delaunayQueue.Count);
            }

            while (this.delaunayQueue.Count > 0)
            {
                HalfEdge he = this.delaunayQueue[0];
                this.delaunayQueue.RemoveAt(0);
                if (he.isType(HalfEdge.AUXILARY))
                {
                    // if (DEBUG) debugView(he, "updateDelaunay: test edge");
                    p1 = he.next.origin;
                    p2 = he.next.next.origin;
                    p3 = he.origin;
                    p4 = he.sibling.next.next.origin;
                    if (this.incircle(p1, p2, p3, p4) > 0 || this.incircle(p3, p4, p1, p2) > 0)
                    {
                        this.flipEdge(he);

                        // if (DEBUG) debugView(he, "updateDelaunay: flipped");
                    }
                }
            }

            // if (TEST) test();
        }

        /// <summary>
        /// The update delaunay all.
        /// </summary>
        public void updateDelaunayAll()
        {
            this.delaunayQueue.AddRange(this.halfEdges);
            this.updateDelaunay();
        }

        /// <summary>
        /// The update interior point.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        public void updateInteriorPoint(Point p)
        {
            //// assert points.contains(p) : error(E_MISSING);
            this.initRemoveConstraints(p);
            this.removeInteriorPoint(p);
            p = this.addInteriorPoint(p);
            p.setType(Point.INTERIOR);
            this.restoreConstraints(p);
            this.updateDelaunay();
        }

        /*    public bool updateBoundaryPoint(BPoint bp, Point2d pNew)
            {
                //// assert points.contains(bp) : error(E_MISSING);

                if (between(bp.prev,bp,pNew)) {
                    return updateBoundaryPointAlong(bp,bp.he,pNew);
                }
                if (between(bp,bp.next,pNew)) {
                    return updateBoundaryPointAlong(bp,bp.next.he,pNew);
                }
                return false;
            }*/

        /// <summary>
        /// The update boundary point along.
        /// </summary>
        /// <param name="bp">
        /// The bp.
        /// </param>
        /// <param name="he">
        /// The he.
        /// </param>
        /// <param name="pNew">
        /// The p new.
        /// </param>
        /// <returns>
        /// The update boundary point along.
        /// </returns>
        private bool updateBoundaryPointAlong(BPoint bp, HalfEdge he, Point2d pNew)
        {
            BPoint bpNew;

            bpNew = new BPoint(pNew);
            this.addBoundaryPoint(bpNew, he);
            this.removeBoundaryPoint(bp, bp.next);
            bp.X = bpNew.X;
            bp.Y = bpNew.Y;
            bp.setType(Point.BOUNDARY);
            bp.prev = bpNew.prev;
            bp.prev.next = bp;
            bp.next = bpNew.next;
            bp.next.prev = bp;
            bp.he = bpNew.he;
            bp.he.origin = bp;

            return true;
        }

        /// <summary>
        /// The update boundary point outside.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        public void updateBoundaryPointOutside(Point p)
        {
            //// assert points.contains(p) : error(E_MISSING);
            int i;
            HalfEdge he;

            if (MESSAGES)
            {
                this.message("Updating boundary point outside existing boundary.");
            }

            this.initRemoveConstraints(p);
            he = p.he.next.next;

            //// assert he.next == p.he : error(E_POLYGON);

            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                if (he.isType(HalfEdge.BOUNDARY))
                {
                    break;
                }

                this.removeEdge(he.sibling);
                /* walk around p counter-clockwise */
                he = he.sibling.next.next;
            }

            //// assert i < halfEdges.Count : E_EXHAUSTED;
            fillGeneralPolygon(p.he);
            this.restoreConstraints(p);
            this.updateDelaunay();

            // if (TEST) test();
        }

        /// <summary>
        /// The update boundary point inside.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <param name="pPrev">
        /// The p prev.
        /// </param>
        /// <param name="pTemp">
        /// The p temp.
        /// </param>
        /// <returns>
        /// The update boundary point inside.
        /// </returns>
        public bool updateBoundaryPointInside(Point p, Point pPrev, Point pTemp)
        {
            //// assert points.contains(p) : error(E_MISSING);
            //// assert points.contains(pPrev) : error(E_MISSING);
            int i;
            Point pNext;
            HalfEdge heTest;

            if (MESSAGES)
            {
                this.message("Updating boundary point inside existing boundary.");
            }

            // locate the previous and next boundary points
            pNext = p.he.next.origin;

            // check for coincidence
            /*if (coincident(p,pTemp)) {
                p.X = pTemp.X;
                p.Y = pTemp.Y;
                return updateBoundaryPointAlong(p,pPrev);
            }*/
            // relocate p to ensure that the quadrilateral (p,pPrev,pTemp,pNext)
            // is a simple polygon
            if (this.intersectProper(pPrev, p, pNext, pTemp))
            {
                if (MESSAGES)
                {
                    this.message("Moving point to ensure simple quadrilateral.");
                }

                Point2d pp = intersection(pPrev, p, pNext, pTemp);
                p.X = pp.X + 0.5 * (pPrev.X - pp.X);
                p.Y = pp.Y + 0.5 * (pPrev.Y - pp.Y);
                this.updateBoundaryPointOutside(p);
            }
            else if (this.intersectProper(pNext, p, pPrev, pTemp))
            {
                if (MESSAGES)
                {
                    this.message("Moving point to ensure simple quadrilateral.");
                }

                Point2d pp = intersection(pNext, p, pPrev, pTemp);
                p.X = pp.X + 0.5 * (pNext.X - pp.X);
                p.Y = pp.Y + 0.5 * (pNext.Y - pp.Y);
                this.updateBoundaryPointOutside(p);
            }

            // insert an interior point where the new boundary point will
            // eventually be
            pTemp = this.addInteriorPoint(pTemp);
            if (MESSAGES)
            {
                this.message("Added temporary interior point.");
            }

            // flip edge pPrev->pNext if it exists
            heTest = pPrev.he;
            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                if (heTest.next.origin == pNext)
                {
                    if (MESSAGES)
                    {
                        this.message("Flipping edge pPrev->pNext.");
                    }

                    this.flipEdge(heTest);
                    break;
                }

                heTest = heTest.next.next;
                if (heTest.isType(HalfEdge.BOUNDARY))
                {
                    break;
                }

                heTest = heTest.sibling;
                if (heTest == pPrev.he)
                {
                    break;
                }
            }

            //// assert i < halfEdges.Count : E_EXHAUSTED;

            // sweep along the new boundary edges to Remove intersecting edges
            this.deleteQueue.Clear();
            this.clearNewBoundaryEdge(pPrev, pTemp);
            this.clearNewBoundaryEdge(pTemp, pNext);

            /* flood fill to Remove points and edges inside the quad */
            if (MESSAGES)
            {
                this.message("Performing flood delete to Remove outside points.");
            }

            var bounds = new Point[4];
            bounds[0] = pPrev;
            bounds[1] = pTemp;
            bounds[2] = pNext;
            bounds[3] = p;
            this.floodDelete(bounds);
            if (MESSAGES)
            {
                this.message("Removing placeholder interior point.");
            }

            this.removePoint(pTemp);
            if (MESSAGES)
            {
                this.message("Relocating boundary point and boundary edges.");
            }

            p.X = pTemp.X;
            p.Y = pTemp.Y;
            fillGeneralPolygon(p.he);
            this.updateDelaunay();
            return true;
        }

        /// <summary>
        /// The validate boundary.
        /// </summary>
        /// <param name="bp">
        /// The bp.
        /// </param>
        /// <returns>
        /// The validate boundary.
        /// </returns>
        private bool validateBoundary(BPoint bp)
        {
            BPoint bpTest;
            bpTest = bp.next;
            for (int j = 3; j < this.nBoundary; j++)
            {
                if (this.intersect(bp.prev, bp, bpTest, bpTest.next))
                {
                    return false;
                }

                bpTest = bpTest.next;
            }

            bpTest = bp.prev;
            for (int j = 3; j < this.nBoundary; j++)
            {
                if (this.intersect(bp, bp.next, bpTest.prev, bpTest))
                {
                    return false;
                }

                bpTest = bpTest.prev;
            }

            return true;
        }

        /// <summary>
        /// The constrain edge.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        /// <returns>
        /// The constrain edge.
        /// </returns>
        private bool constrainEdge(HalfEdge he)
        {
            //// assert halfEdges.contains(he) : error(E_MISSING);
            if (MESSAGES)
            {
                this.message("Constraining edge %d.", this.halfEdges.IndexOf(he));
            }

            if (he.isType(HalfEdge.BOUNDARY))
            {
                if (MESSAGES)
                {
                    this.message("Ignoring boundary edge constraint.");
                }

                return false;
            }

            he.constrain();
            he.sibling.constrain();
            return true;
        }

        /// <summary>
        /// The constrain all edges.
        /// </summary>
        public void constrainAllEdges()
        {
            foreach (HalfEdge he in this.halfEdges)
            {
                if (!he.isType(HalfEdge.BOUNDARY))
                {
                    he.constrain();
                }
            }
        }

        /// <summary>
        /// The init remove constraints.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        public void initRemoveConstraints(Point p)
        {
            //// assert points.contains(p) : error(E_MISSING);
            this.removeConstraintPeg = p;
            this.removedConstraints.Clear();
        }

        /// <summary>
        /// The restore constraints.
        /// </summary>
        /// <param name="p">
        /// The p.
        /// </param>
        public void restoreConstraints(Point p)
        {
            //// assert points.contains(p) : error(E_MISSING);
            //// assert p == removeConstraintPeg : error("Race condition!");
            this.removedConstraints.Remove(p);
            foreach (Point p0 in this.removedConstraints)
            {
                if (p0.isType(Point.DELETED))
                {
                    continue;
                }

                this.addConstraint(p, p0);
            }

            this.removeConstraintPeg = null;
        }

        /// <summary>
        /// The flip edge.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        /// <returns>
        /// The flip edge.
        /// </returns>
        private bool flipEdge(HalfEdge he)
        {
            //// assert halfEdges.contains(he) : error(E_MISSING);
            //// assert ! he.isType(HalfEdge.BOUNDARY) : error(E_TYPE);
            HalfEdge he1, he2, he3, he4;

            if (MESSAGES)
            {
                this.message("Flipping edge %d.", this.halfEdges.IndexOf(he));
            }

            // locate halfedges
            he1 = he.next;
            he2 = he1.next;
            he3 = he.sibling.next;
            he4 = he3.next;

            // flip the origins
            he.origin = he2.origin;
            he.sibling.origin = he4.origin;

            // update point->halfedge pointers
            this.updateHalfEdge(he3);
            this.updateHalfEdge(he1);

            // link halfedges
            he1.next = he;
            he.next = he4;
            he4.next = he1;
            he3.next = he.sibling;
            he.sibling.next = he2;
            he2.next = he3;

            // Add halfedges to the delaunay test
            this.delaunayQueue.Add(he1);
            this.delaunayQueue.Add(he2);
            this.delaunayQueue.Add(he3);
            this.delaunayQueue.Add(he4);
            return true;
        }

        /******************************************************************************/
        /* Locator methods
        /******************************************************************************/

        /**
         * This brute force approach works for *any* non-intersecting boundary,
         * concave or convex. If the boundary is guaranteed to be a convex, a
         * smarter face-walking algorithm could be used.
         */

        /// <summary>
        /// The find face brute force.
        /// </summary>
        /// <param name="heStart">
        /// The he start.
        /// </param>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <returns>
        /// </returns>
        public FaceWalk findFaceBruteForce(HalfEdge heStart, Point p)
        {
            var ccw = new double[3];
            HalfEdge he1, he2;

            this.clearFlags(HalfEdge.FLAG_ALGORITHM);
            foreach (HalfEdge he0 in this.halfEdges)
            {
                if (he0.isFlagged(HalfEdge.FLAG_ALGORITHM))
                {
                    continue;
                }

                he1 = he0.next;
                he2 = he1.next;

                //// assert he2.next == he0 : error("Found non-face!");
                he0.flag(HalfEdge.FLAG_ALGORITHM);
                he1.flag(HalfEdge.FLAG_ALGORITHM);
                he2.flag(HalfEdge.FLAG_ALGORITHM);
                ccw[0] = this.orient(he0.origin, he1.origin, p);
                if (ccw[0] < 0)
                {
                    continue;
                }

                ccw[1] = this.orient(he1.origin, he2.origin, p);
                if (ccw[1] < 0)
                {
                    continue;
                }

                ccw[2] = this.orient(he2.origin, he0.origin, p);
                if (ccw[2] < 0)
                {
                    continue;
                }

                if (ccw[0] == 0)
                {
                    return new FaceWalk(he0, FaceWalk.COINCIDENT);
                }

                if (ccw[1] == 0)
                {
                    return new FaceWalk(he1, FaceWalk.COINCIDENT);
                }

                if (ccw[2] == 0)
                {
                    return new FaceWalk(he2, FaceWalk.COINCIDENT);
                }

                return new FaceWalk(he0, FaceWalk.CLOCKWISE);
            }

            if (DEBUG)
            {
                this.error("Exhausted faces!");
            }

            return null;
        }

        /**
         * A slightly smarter face walk routine that resorts to brute force
         * only when it gets confused by an concave boundary.
         */

        /// <summary>
        /// The find face.
        /// </summary>
        /// <param name="heStart">
        /// The he start.
        /// </param>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <returns>
        /// </returns>
        /// <exception cref="Exception">
        /// </exception>
        public FaceWalk findFace(HalfEdge heStart, Point p)
        {
            int i;
            var ccw = new double[3];
            HalfEdge he0, he1, he2;
            var queue = new List<HalfEdge>();

            this.clearFlags(HalfEdge.FLAG_ALGORITHM);
            queue.Add(heStart);
            queue.AddRange(this.halfEdges);

            // he0 = queue.pop();
            he0 = queue[0];
            queue.RemoveAt(0);
            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                if (he0 == null)
                {
                    throw new Exception();
                }

                if (he0.isFlagged(HalfEdge.FLAG_ALGORITHM))
                {
                    he0 = queue[0];
                    queue.RemoveAt(0);
                    continue;
                }

                he1 = he0.next;
                he2 = he1.next;

                //// assert he2.next == he0 : error("Found non-face!");
                he0.flag(HalfEdge.FLAG_ALGORITHM);
                he1.flag(HalfEdge.FLAG_ALGORITHM);
                he2.flag(HalfEdge.FLAG_ALGORITHM);
                ccw[0] = this.orient(he0.origin, he1.origin, p);
                if (ccw[0] < 0)
                {
                    he0 = he0.sibling;
                    continue;
                }

                ccw[1] = this.orient(he1.origin, he2.origin, p);
                if (ccw[1] < 0)
                {
                    he0 = he1.sibling;
                    continue;
                }

                ccw[2] = this.orient(he2.origin, he0.origin, p);
                if (ccw[2] < 0)
                {
                    he0 = he2.sibling;
                    continue;
                }

                if (ccw[0] == 0)
                {
                    return new FaceWalk(he0, FaceWalk.COINCIDENT);
                }

                if (ccw[1] == 0)
                {
                    return new FaceWalk(he1, FaceWalk.COINCIDENT);
                }

                if (ccw[2] == 0)
                {
                    return new FaceWalk(he2, FaceWalk.COINCIDENT);
                }

                return new FaceWalk(he0, FaceWalk.CLOCKWISE);
            }

            //// assert i < halfEdges.Count : error(E_EXHAUSTED);
            if (DEBUG)
            {
                this.error("Exhausted faces!");
            }

            return null;
        }

        /// <summary>
        /// The find previous.
        /// </summary>
        /// <param name="he">
        /// The he.
        /// </param>
        /// <returns>
        /// </returns>
        protected HalfEdge findPrevious(HalfEdge he)
        {
            //// assert halfEdges.contains(he) : error(E_MISSING,he);
            int i;
            HalfEdge heSearch;

            heSearch = he.next;
            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                // if (DEBUG) debugView(heSearch, "findPrevious: searching forward");
                if (heSearch.next == he)
                {
                    break;
                }

                heSearch = heSearch.next;
            }

            //// assert i < halfEdges.Count : error(E_EXHAUSTED,he);

            //// assert halfEdges.contains(heSearch) : error(E_MISSING,heSearch);
            return heSearch;
        }

        // robust with non-triangular regions
        /// <summary>
        /// The start face walk.
        /// </summary>
        /// <param name="pStart">
        /// The p start.
        /// </param>
        /// <param name="pEnd">
        /// The p end.
        /// </param>
        /// <returns>
        /// </returns>
        private FaceWalk startFaceWalk(Point pStart, Point pEnd)
        {
            //// assert points.contains(pStart) : error(E_MISSING);
            //// assert points.contains(pEnd) : error(E_MISSING);
            //// assert pStart != pEnd : error(E_IDENTICAL,pStart,pEnd);
            //// assert ! coincident(pStart,pEnd) : error(E_COINCIDENT,pStart,pEnd);
            int i;
            double ccwTrailing, ccwLeading;
            Point pTrailing, pLeading;
            HalfEdge he;
            HalfEdge hePrev = null;

            he = pStart.he;
            pTrailing = he.next.origin;
            ccwTrailing = this.orient(pStart, pEnd, pTrailing);

            // special case for boundary starting points
            if (pStart.isType(Point.BOUNDARY))
            {
                // check for coincidence with boundary
                if (pTrailing == pEnd)
                {
                    return new FaceWalk(he, FaceWalk.COINCIDENT);
                }

                // check whether pStart->pEnd is within epsilon of the boundary
                //// assert ccwTrailing <= 0 :
                // error("Point lies outside boundary!", pEnd);
                // it is safe to assume that he is clockwise because of the above
                // //// assertion
                if (ccwTrailing == 0)
                {
                    if (this.betweenProper(pStart, pEnd, pTrailing) || this.betweenProper(pStart, pTrailing, pEnd))
                    {
                        return new FaceWalk(he, FaceWalk.CLOCKWISE);
                    }
                }
            }

            for (i = 0; i <= this.halfEdges.Count; i++)
            {
                // this face may be a polygon, so search forward to find
                // the second test edge
                hePrev = this.findPrevious(he);
                pTrailing = he.next.origin;

                // if (DEBUG) debugView(pTrailing, "startFaceWalk: trailing point");
                pLeading = hePrev.origin;

                // if (DEBUG) debugView(pLeading, "startFaceWalk: leading point");
                // check for coincidence of either star edge
                if (pTrailing == pEnd)
                {
                    return new FaceWalk(he, FaceWalk.COINCIDENT);
                }
                else if (pLeading == pEnd)
                {
                    return new FaceWalk(hePrev, FaceWalk.COINCIDENT);
                }

                //// assert ! coincident(pTrailing,pEnd) :
                // error(E_COINCIDENT, pTrailing, pEnd);
                //// assert ! coincident(pLeading,pEnd) :
                // error(E_COINCIDENT, pLeading, pEnd);
                // check if the leading point is counter-clockwise/collinear and the
                // trailing point clockwise of pStart->pEnd
                ccwLeading = this.orient(pStart, pEnd, pLeading);
                if (ccwLeading >= 0 && ccwTrailing < 0)
                {
                    return new FaceWalk(he, FaceWalk.CLOCKWISE);
                }

                ccwTrailing = ccwLeading;
                he = hePrev.sibling;
            }

            //// assert i < halfEdges.Count : E_EXHAUSTED;

            // return a failed walk at this point
            return new FaceWalk(null, FaceWalk.FAILED);
        }

        /******************************************************************************/
        /* Drawing methods
        /******************************************************************************/

        /*
        public void drawDebug(Graphics2D g2)
    {
            
        if (DEBUG) {
            g2.setColor(highlightColor);
            foreach (Object obj in debugObjects) {
                if (obj instanceof Point) {
                    drawPoint(g2,(Point)obj);
                } else if (obj instanceof HalfEdge) {
                    drawHalfEdge(g2,(HalfEdge)obj);
                    drawLine(g2,(HalfEdge)obj);
                } else if (obj instanceof Edge) {
                    Edge e = (Edge)obj;
                    drawHalfEdge(g2,e.p2,e.p1);
                    drawLine(g2,e.p1,e.p2);
                }
            }
        }
             
    }
        */
        /******************************************************************************/
        /* Computational geometry methods and predicates
        /******************************************************************************/

        // compute the normalized projection of ac onto ab
        /// <summary>
        /// The proj norm.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <param name="c">
        /// The c.
        /// </param>
        /// <returns>
        /// The proj norm.
        /// </returns>
        private static double projNorm(Point2d a, Point2d b, Point2d c)
        {
            double x1, x2, y1, y2;
            x1 = b.X - a.X;
            x2 = c.X - a.X;
            y1 = b.Y - a.Y;
            y2 = c.Y - a.Y;
            return (x1 * x2 + y1 * y2) / (x1 * x1 + y1 * y1);
        }

        // compute the magnitude of the cross product of ab and ac
        /// <summary>
        /// The cross.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <param name="c">
        /// The c.
        /// </param>
        /// <returns>
        /// The cross.
        /// </returns>
        private static double cross(Point2d a, Point2d b, Point2d c)
        {
            double x1, x2, y1, y2;
            x1 = b.X - a.X;
            x2 = c.X - a.X;
            y1 = b.Y - a.Y;
            y2 = c.Y - a.Y;
            return x1 * y2 - y1 * x2;
        }

        // compute the squared perpendicular distance of c onto ab
        /// <summary>
        /// The perp dist sq.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <param name="c">
        /// The c.
        /// </param>
        /// <returns>
        /// The perp dist sq.
        /// </returns>
        private static double perpDistSq(Point2d a, Point2d b, Point2d c)
        {
            double x1, x2, y1, y2, cross, lenSq;
            x1 = b.X - a.X;
            x2 = c.X - a.X;
            y1 = b.Y - a.Y;
            y2 = c.Y - a.Y;
            cross = x1 * y2 - y1 * x2;
            lenSq = cross * cross;
            lenSq /= x1 * x1 + y1 * y1;
            return lenSq;
        }

        /// <summary>
        /// The coincident.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <returns>
        /// The coincident.
        /// </returns>
        public bool coincident(Point2d a, Point2d b)
        {
            if (FindSquareDist(a as Point, b as Point) < this.epsilon)
            {
                return true;
            }

            return false;
        }

        // project the point p orthogonally onto the segment p1->p2
        // and return the scaled, parameterized position in [0,1]
        /// <summary>
        /// The projection.
        /// </summary>
        /// <param name="p1">
        /// The p 1.
        /// </param>
        /// <param name="p2">
        /// The p 2.
        /// </param>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <returns>
        /// The projection.
        /// </returns>
        public static double projection(Point2d p1, Point2d p2, Point2d p)
        {
            double ax, ay, bx, by;
            ax = p.X - p1.X;
            ay = p.Y - p1.Y;
            bx = p2.X - p1.X;
            by = p2.Y - p1.Y;
            return (ax * bx + ay * by) / (bx * bx + by * by);
        }

        // returns the intersection point of segments ab and cd
        // as a point on ab
        /// <summary>
        /// The intersection.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <param name="c">
        /// The c.
        /// </param>
        /// <param name="d">
        /// The d.
        /// </param>
        /// <returns>
        /// </returns>
        public static Point2d intersection(Point2d a, Point2d b, Point2d c, Point2d d)
        {
            double t, l1, l2;
            double cdx, cdy;
            Point2d p;

            cdx = c.X - d.X;
            cdy = c.Y - d.Y;

            // distance from a to cd
            l1 = Math.Abs((a.X - d.X) * cdy - (a.Y - d.Y) * cdx);

            // distance from b to cd
            l2 = Math.Abs((b.X - d.X) * cdy - (b.Y - d.Y) * cdx);

            // need to handle case where l1+l2 = 0
            // if this method could be called on parallel segments
            // that overlap
            if (l1 + l2 == 0)
            {
                // System.err.println(                        "Intersection called on parallel overlapping segments!");
            }

            t = l1 / (l1 + l2);
            p = new Point2d((1 - t) * a.X + t * b.X, (1 - t) * a.Y + t * b.Y);
            return p;
        }

        // from O'Rourke's Computational Geometry in C
        /// <summary>
        /// The intersect.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <param name="c">
        /// The c.
        /// </param>
        /// <param name="d">
        /// The d.
        /// </param>
        /// <returns>
        /// The intersect.
        /// </returns>
        public bool intersect(Point2d a, Point2d b, Point2d c, Point2d d)
        {
            if (this.intersectProper(a, b, c, d))
            {
                return true;
            }
            else if (this.between(a, b, c) || this.between(a, b, d) || this.between(c, d, a) || this.between(c, d, b))
            {
                return true;
            }
            else
            {
                return false;
            }
        }

        // from O'Rourke's Computational Geometry in C
        /// <summary>
        /// The intersect proper.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <param name="c">
        /// The c.
        /// </param>
        /// <param name="d">
        /// The d.
        /// </param>
        /// <returns>
        /// The intersect proper.
        /// </returns>
        public bool intersectProper(Point2d a, Point2d b, Point2d c, Point2d d)
        {
            /* Eliminate improper cases. */
            if (this.orient(a, b, c) == 0 || this.orient(a, b, d) == 0 || this.orient(c, d, a) == 0 ||
                this.orient(c, d, b) == 0)
            {
                return false;
            }
            else if (this.orient(a, b, c) * this.orient(a, b, d) > 0 || this.orient(c, d, a) * this.orient(c, d, b) > 0)
            {
                return false;
            }
            else
            {
                return true;
            }
        }

        /*
         * Tests whether c is within the epsilon tubular neighborhood
         * around segment ab.
         */

        /// <summary>
        /// The between.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <param name="c">
        /// The c.
        /// </param>
        /// <returns>
        /// The between.
        /// </returns>
        public bool between(Point2d a, Point2d b, Point2d c)
        {
            /* check the epsilon neighborhood at the endpoints */
            if (this.coincident(a, c))
            {
                return true;
            }
            else if (this.coincident(b, c))
            {
                return true;
            }
            else
            {
                /* check the epsilon neighborhood along the segment */
                if (perpDistSq(a, b, c) < this.epsilon)
                {
                    double d = projNorm(a, b, c);
                    if (0 < d && d < 1)
                    {
                        return true;
                    }
                }
            }

            return false;
        }

        /*
         * Tests whether c is within the epsilon tubular neighborhood
         * around segment ab, but excludes the epsilon neighborhoods
         * around a and b.
         */

        /// <summary>
        /// The between proper.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <param name="c">
        /// The c.
        /// </param>
        /// <returns>
        /// The between proper.
        /// </returns>
        public bool betweenProper(Point2d a, Point2d b, Point2d c)
        {
            /* reject the epsilon neighborhood at the endpoints */
            if (this.coincident(a, c))
            {
                return false;
            }
            else if (this.coincident(b, c))
            {
                return false;
            }
            else
            {
                /* check the epsilon neighborhood along the segment */
                if (perpDistSq(a, b, c) < this.epsilon)
                {
                    double d = projNorm(a, b, c);
                    if (0 < d && d < 1)
                    {
                        return true;
                    }
                }
            }

            return false;
        }

        /*
         * Orientation predicate that returns positive if ray a->b must
         * turn counter-clockwise to intersect vertex c.
         *
         * Adapted from Jonathan Shewchuk's predicates.c orient2d() routine.
         * This implementation uses only a single-level floating-point filter,
         * which is a simplified and less efficient version of his multi-level
         * adaptive method.
         */

        /// <summary>
        /// The orient.
        /// </summary>
        /// <param name="pa">
        /// The pa.
        /// </param>
        /// <param name="pb">
        /// The pb.
        /// </param>
        /// <param name="pc">
        /// The pc.
        /// </param>
        /// <returns>
        /// The orient.
        /// </returns>
        public double orient(Point2d pa, Point2d pb, Point2d pc)
        {
            double detleft, detright, det;
            double detsum, errbound;

            detleft = (pa.X - pc.X) * (pb.Y - pc.Y);
            detright = (pa.Y - pc.Y) * (pb.X - pc.X);
            det = detleft - detright;

            if (detleft > 0.0)
            {
                if (detright <= 0.0)
                {
                    return det;
                }
                else
                {
                    detsum = detleft + detright;
                }
            }
            else if (detleft < 0.0)
            {
                if (detright >= 0.0)
                {
                    return det;
                }
                else
                {
                    detsum = -detleft - detright;
                }
            }
            else
            {
                return det;
            }

            errbound = this.orientErrorBound * detsum;
            if ((det >= errbound) || (-det >= errbound))
            {
                return det;
            }

            if (MESSAGES)
            {
                this.message("orient = %g\n", det);
            }

            /* resort to exact arithmetic */
            return this.orientExact(pa, pb, pc);
        }

        /*
         * Orientation predicate that returns positive if ray a->b must
         * turn counter-clockwise to intersect vertex c.
         *
         * Performs an exact calculation using the Apfloat arbitrary-precision
         * library.
         */

        /// <summary>
        /// The orient exact.
        /// </summary>
        /// <param name="pa">
        /// The pa.
        /// </param>
        /// <param name="pb">
        /// The pb.
        /// </param>
        /// <param name="pc">
        /// The pc.
        /// </param>
        /// <returns>
        /// The orient exact.
        /// </returns>
        public double orientExact(Point2d pa, Point2d pb, Point2d pc)
        {
            float ax, ay, bx, by, cx, cy;
            float acx, bcx, acy, bcy;
            float detleft, detright, det;

            det = 0;
            {
                ax = (float)pa.X;
                ay = (float)pa.Y;
                bx = (float)pb.X;
                by = (float)pb.Y;
                cx = (float)pc.X;
                cy = (float)pc.Y;

                acx = ax - cx;
                bcx = bx - cx;
                acy = ay - cy;
                bcy = by - cy;

                detleft = acx * bcy;
                detright = acy * bcx;

                det = detleft - detright;
            }

            // if (MESSAGES) message("orientExact = %g\n", det.doubleValue());
            return det.CompareTo(0);
        }

        /*
         * Non-robust, deprecated orientation predicate.
         */

        /// <summary>
        /// The orient non robust.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <param name="c">
        /// The c.
        /// </param>
        /// <returns>
        /// The orient non robust.
        /// </returns>
        public double orientNonRobust(Point2d a, Point2d b, Point2d c)
        {
            if (perpDistSq(a, b, c) < this.epsilon)
            {
                return 0;
            }

            return cross(a, b, c);
        }

        /*
         * Incircle predicate that returns positive if point pd lies inside
         * the circumcircle through pa, pb, and pc.
         *
         * Adapted from Jonathan Shewchuk's predicates.c incircle() routine.
         * This implementation uses only a single-level floating-point filter,
         * which is a simplified and less efficient version of his multi-level
         * adaptive method.
         */

        /// <summary>
        /// The incircle.
        /// </summary>
        /// <param name="pa">
        /// The pa.
        /// </param>
        /// <param name="pb">
        /// The pb.
        /// </param>
        /// <param name="pc">
        /// The pc.
        /// </param>
        /// <param name="pd">
        /// The pd.
        /// </param>
        /// <returns>
        /// The incircle.
        /// </returns>
        public double incircle(Point2d pa, Point2d pb, Point2d pc, Point2d pd)
        {
            double adx, ady, bdx, bdy, cdx, cdy;
            double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
            double alift, blift, clift;
            double det;
            double permanent, errbound;

            adx = pa.X - pd.X;
            bdx = pb.X - pd.X;
            cdx = pc.X - pd.X;
            ady = pa.Y - pd.Y;
            bdy = pb.Y - pd.Y;
            cdy = pc.Y - pd.Y;

            bdxcdy = bdx * cdy;
            cdxbdy = cdx * bdy;
            alift = adx * adx + ady * ady;

            cdxady = cdx * ady;
            adxcdy = adx * cdy;
            blift = bdx * bdx + bdy * bdy;

            adxbdy = adx * bdy;
            bdxady = bdx * ady;
            clift = cdx * cdx + cdy * cdy;

            det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady);

            if (bdxcdy < 0)
            {
                bdxcdy = -bdxcdy;
            }

            if (cdxbdy < 0)
            {
                cdxbdy = -cdxbdy;
            }

            if (cdxady < 0)
            {
                cdxady = -cdxady;
            }

            if (adxcdy < 0)
            {
                adxcdy = -adxcdy;
            }

            if (adxbdy < 0)
            {
                adxbdy = -adxbdy;
            }

            if (bdxady < 0)
            {
                bdxady = -bdxady;
            }

            permanent = (bdxcdy + cdxbdy) * alift + (cdxady + adxcdy) * blift + (adxbdy + bdxady) * clift;
            errbound = this.incircleErrorBound * permanent;
            if ((det > errbound) || (-det > errbound))
            {
                return det;
            }

            return this.incircleExact(pa, pb, pc, pd);
        }

        /*
         * Incircle predicate that returns positive if point pd lies inside
         * the circumcircle through pa, pb, and pc.
         *
         * Performs an exact calculation using the Apfloat arbitrary-precision
         * library.
         */

        /// <summary>
        /// The incircle exact.
        /// </summary>
        /// <param name="pa">
        /// The pa.
        /// </param>
        /// <param name="pb">
        /// The pb.
        /// </param>
        /// <param name="pc">
        /// The pc.
        /// </param>
        /// <param name="pd">
        /// The pd.
        /// </param>
        /// <returns>
        /// The incircle exact.
        /// </returns>
        public double incircleExact(Point2d pa, Point2d pb, Point2d pc, Point2d pd)
        {
            float ax, ay, bx, by, cx, cy, dx, dy;
            float adx, ady, bdx, bdy, cdx, cdy;
            float bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
            float alift, blift, clift;
            float det;

            det = 0;
            {
                ax = (float)pa.X;
                ay = (float)pa.Y;
                bx = (float)pb.X;
                by = (float)pb.Y;
                cx = (float)pc.X;
                cy = (float)pc.Y;
                dx = (float)pd.X;
                dy = (float)pd.Y;

                dx = -dx;
                dy = -dy;

                adx = ax + dx;
                bdx = bx + dx;
                cdx = cx + dx;
                ady = ay + dy;
                bdy = by + dy;
                cdy = cy + dy;

                bdxcdy = bdx * cdy;
                cdxbdy = cdx * bdy;

                cdxady = cdx * ady;
                adxcdy = adx * cdy;

                adxbdy = adx * bdy;
                bdxady = bdx * ady;

                adx = adx * adx;
                ady = ady * ady;
                alift = adx + ady;

                bdx = bdx * bdx;
                bdy = bdy * bdy;
                blift = bdx + bdy;

                cdx = cdx * cdx;
                cdy = cdy * cdy;
                clift = cdx + cdy;

                alift = alift * (bdxcdy - cdxbdy);
                blift = blift * (cdxady - adxcdy);
                clift = clift * (adxbdy - bdxady);

                det = alift + blift + clift;
            }

            // if (MESSAGES) message("incircleExact = %g\n", det.doubleValue());
#warning сравнения не понятно
            return det.CompareTo(0);
        }

        /*
         * Non-robust, deprecated incircle predicate.
         */

        /// <summary>
        /// The incircle non robust.
        /// </summary>
        /// <param name="a">
        /// The a.
        /// </param>
        /// <param name="b">
        /// The b.
        /// </param>
        /// <param name="c">
        /// The c.
        /// </param>
        /// <param name="d">
        /// The d.
        /// </param>
        /// <returns>
        /// The incircle non robust.
        /// </returns>
        public static int incircleNonRobust(Point2d a, Point2d b, Point2d c, Point2d d)
        {
            double adx, ady, bdx, bdy, cdx, cdy;
            double abdet, bcdet, cadet;
            double alift, blift, clift;

            adx = a.X - d.X;
            ady = a.Y - d.Y;
            bdx = b.X - d.X;
            bdy = b.Y - d.Y;
            cdx = c.X - d.X;
            cdy = c.Y - d.Y;

            abdet = adx * bdy - bdx * ady;
            bcdet = bdx * cdy - cdx * bdy;
            cadet = cdx * ady - adx * cdy;
            alift = adx * adx + ady * ady;
            blift = bdx * bdx + bdy * bdy;
            clift = cdx * cdx + cdy * cdy;

            return Math.Sign(alift * bcdet + blift * cadet + clift * abdet);
        }

        // adapted from java.awt.geom.Line2d.ptSegDistSq()
        /// <summary>
        /// The edge distance sq.
        /// </summary>
        /// <param name="p1">
        /// The p 1.
        /// </param>
        /// <param name="p2">
        /// The p 2.
        /// </param>
        /// <param name="p">
        /// The p.
        /// </param>
        /// <returns>
        /// The edge distance sq.
        /// </returns>
        public static double edgeDistanceSq(Point2d p1, Point2d p2, Point2d p)
        {
            double x2, y2, px, py;

            // Adjust vectors relative to x1,y1
            // x2,y2 becomes relative vector from x1,y1 to end of segment
            x2 = p2.X - p1.X;
            y2 = p2.Y - p1.Y;

            // px,py becomes relative vector from x1,y1 to test point
            px = p.X - p1.X;
            py = p.Y - p1.Y;
            double dotprod = px * x2 + py * y2;
            double projlenSq;
            if (dotprod <= 0.0)
            {
                // px,py is on the side of x1,y1 away from x2,y2
                // distance to segment is length of px,py vector
                // "length of its (clipped) projection" is now 0.0
                projlenSq = 0.0;
            }
            else
            {
                // switch to backwards vectors relative to x2,y2
                // x2,y2 are already the negative of x1,y1=>x2,y2
                // to get px,py to be the negative of px,py=>x2,y2
                // the dot product of two negated vectors is the same
                // as the dot product of the two normal vectors
                px = x2 - px;
                py = y2 - py;
                dotprod = px * x2 + py * y2;
                if (dotprod <= 0.0)
                {
                    // px,py is on the side of x2,y2 away from x1,y1
                    // distance to segment is length of (backwards) px,py vector
                    // "length of its (clipped) projection" is now 0.0
                    projlenSq = 0.0;
                }
                else
                {
                    // px,py is between x1,y1 and x2,y2
                    // dotprod is the length of the px,py vector
                    // projected on the x2,y2=>x1,y1 vector times the
                    // length of the x2,y2=>x1,y1 vector
                    projlenSq = dotprod * dotprod / (x2 * x2 + y2 * y2);
                }
            }

            // Distance to line is now the length of the relative point
            // vector minus the length of its projection onto the line
            // (which is zero if the projection falls outside the range
            // of the line segment).
            double lenSq = px * px + py * py - projlenSq;
            if (lenSq < 0)
            {
                lenSq = 0;
            }

            return lenSq;
        }

        /******************************************************************************/
        /* Debug, test, and output methods
        /******************************************************************************/
        /*
        public void test()
        {
            message("Testing.");
            testing = true;
            testPoints();
            testFaces();
            testHalfEdges();
            testHalfEdgePointers();
            testing = false;
        }

        private void testPoints()
        {
            for (int i = 0; i < points.Count; i++)
            {
                Point p1 = points[i);
                for (int j = 0; j < points.Count; j++)
                {
                    Point p2 = points[j);
                    if (p1 != p2 && coincident(p1, p2))
                    {
                        System.err.printf("%s: Coincident points: %d, %d!\n",
                                name, i, j);
                    }
                }
                for (int j=0; j<halfEdges.Count; j++)  {
                    HalfEdge he = halfEdges[j);
                    if (between(he.origin,he.next.origin,p1) &&
                        he.origin != p1 &&
                        he.next.origin != p1)
                    {
                        System.err.printf("%s: Point %d overlaps halfedge %d!\n",
                                name,i,j);
                        if (DEBUG) debugView(p1,"Point overlaps halfedge");
                        if (DEBUG) debugView(he,"Point overlaps halfedge");
                    }
                }
            }
        }*/
        /*
        private void testFaces()
    {
        int i;
        HalfEdge heTest;
        LinkedList<HalfEdge> used = new LinkedList();
        LinkedList<HalfEdge> face = new LinkedList();
        
        foreach (HalfEdge he in halfEdges) {
            if (used.contains(he)) continue;
            face.Clear();
            heTest = he;
            for (i=0; i<=halfEdges.Count; i++) {
                used.Add(heTest);
                face.Add(heTest);
                heTest = heTest.next;
                if (heTest == he) break;
            } //// assert i < halfEdges.Count : error(E_EXHAUSTED);
            if (face.Count > 3) {
                System.out.printf("%s: polygon: (",name);
                for (HalfEdge he0 : face) {
                        System.out.printf(" %d ",points.IndexOf(he0.origin));
                    }
                    System.out.println(")");
                }
            }
        }
        
        private void testHalfEdges()
        {
            for (int i = 0; i < halfEdges.Count; i++)
            {
                HalfEdge he = halfEdges[i);
                if (points.contains(he.origin) == false)
                {
                    System.err.printf(
                            "%s: Missing origin for halfedge %d!\n",
                            name, i);
                }
                if (halfEdges.contains(he.next) == false)
                {
                    System.err.printf(
                            "%s: Missing next pointer for halfedge %d!\n",
                            name, i);
                }
                if (coincident(he.origin, he.next.origin))
                {
                    System.err.printf(
                            "%s: Coincident endpoints for halfedge %d!\n",
                            name, i);
                }
                if (!he.isType(HalfEdge.BOUNDARY))
                {
                    if (halfEdges.contains(he.sibling) == false)
                    {
                        System.err.printf(
                                "%s: Missing sibling for halfedge %d!\n",
                                name, i);
                    }
                    if (he.sibling.sibling != he)
                    {
                        System.err.printf(
                                "%s: Mismatched sibling for halfedge %d!\n",
                                name, i);
                    }
                    if (he.next.origin != he.sibling.origin)
                    {
                        System.err.printf(
                                "%s: Unaligned next/sibling for halfedge %d\n",
                                name, i);
                    }
                }
            }
        }

        private void testHalfEdgePointers()
        {
            for (int i = 0; i < points.Count; i++)
            {
                Point p = points[i);
                if (halfEdges.contains(p.he) == false)
                {
                    System.err.printf(
                            "%s: Missing halfedge for point %d!\n",
                            name, i);
                }
                if (p.he.origin != p)
                {
                    System.err.printf(
                            "%s: Mismatched halfedge for point %d!\n",
                            name, i);
                }
            }
        }
        */
        /*
      public void listPoints()
      {
          message("### POINT LIST ###");
          for (int i = 0; i < points.Count; i++)
          {
              Point p = points[i);
              if (i % 20 == 0)
              {
                  message("     ID | Halfedge |    Pair | Type");
              }
              message("%7d |  %7d | %7d |    %1d",
                          i,
                          halfEdges.IndexOf(p.he),
                          points.IndexOf(p.pair),
                          p.type);
          }
      }*/
        /*
        public void listHalfEdges()
        {
            message("### HALFEDGE LIST ###");
            for (int i = 0; i < halfEdges.Count; i++)
            {
                HalfEdge he = halfEdges[i);
                if (i % 20 == 0)
                {
                    message("     ID ->    Next |  Origin ->    " +
                            "Next | Sibling | Type");
                }
                message("%7d -> %7d | %7d -> %7d | %7d | %1d",
                        i,
                        halfEdges.IndexOf(he.next),
                        points.IndexOf(he.origin),
                        points.IndexOf(he.next.origin),
                        halfEdges.IndexOf(he.sibling),
                        he.type);
            }
        }*/
        /*
        public void debugView()
        {
            debugView("(no message)");
        }
    
        public void debugView(String s)
        {
            if (toggleDebugVisual) {
                //debugFrame.repaint();
                debugPanel.repaint();
            }
            if (toggleDebugInteractive
                && !toggleDebugSuspend )
            {
                // pause for user response
                JOptionPane pane = new JOptionPane(
                        s + "\nContinue?\n",
                        JOptionPane.QUESTION_MESSAGE,
                        JOptionPane.YES_NO_CANCEL_OPTION);
                JDialog dialog = pane.createDialog(debugFrame,"Paused");
                dialog.setLocation(0,0);
                dialog.setVisible(true);
                Integer selectedValue = (Integer)pane.getValue();
                if (selectedValue == JOptionPane.CANCEL_OPTION) {
                    toggleDebugSuspend = true;
                } else if (selectedValue == JOptionPane.NO_OPTION) {
                    System.exit(0);
                }
            }
        }
         */
        /*
        public void debugView(Point p, String s)
        {
            debugObjects.Add(p);
            debugView(s);
            debugObjects.Remove(p);
        }
    
        public void debugView(HalfEdge he, String s)
        {
            debugObjects.Add(he);
            debugView(s);
            debugObjects.Remove(he);
        }
        */
        /*
        public void debugView(String s, Object ... args)
        {
            for (Object obj : args) {
                debugObjects.Add(obj);
            }
            debugView(s);
            for (Object obj : args) {
                debugObjects.Remove(obj);
            }
        }
        */

        /// <summary>
        /// The dv start.
        /// </summary>
        public void dvStart()
        {
            this.toggleDebugInteractive = this.toggleDebugVisual = true;
        }

        /// <summary>
        /// The dv end.
        /// </summary>
        public void dvEnd()
        {
            this.toggleDebugInteractive = this.toggleDebugVisual = false;
        }

        /*
        public void debugPause(String s, Object ... args)
        {
            dvStart();
            for (Object obj : args) {
                debugView(s,obj);
            }
            dvEnd();
        }
        */

        /// <summary>
        /// The message.
        /// </summary>
        /// <param name="s">
        /// The s.
        /// </param>
        public void message(string s)
        {
            Debug.Write(": ");
            Debug.Write(s);
            Debug.Write("\n");
        }

        /// <summary>
        /// The message.
        /// </summary>
        /// <param name="s">
        /// The s.
        /// </param>
        /// <param name="args">
        /// The args.
        /// </param>
        public void message(string s, object args)
        {
            this.message(String.Format(s));
        }

        /// <summary>
        /// The message.
        /// </summary>
        /// <param name="s">
        /// The s.
        /// </param>
        /// <param name="args1">
        /// The args 1.
        /// </param>
        /// <param name="args2">
        /// The args 2.
        /// </param>
        public void message(string s, object args1, object args2)
        {
            this.message(String.Format(s));
        }

        // used in //// assert statements to dump halfedge and point lists

        /// <summary>
        /// The error.
        /// </summary>
        /// <param name="s">
        /// The s.
        /// </param>
        /// <returns>
        /// The error.
        /// </returns>
        protected string error(string s)
        {
            // listHalfEdges();
            // listPoints();
            // if (!testing) test();
            return s;
        }

        #region  Drawing methods

        /*
        public void drawPoints(Graphics2D g2)
        {
            foreach (Point p in points)
            {
                g2.setColor(p.getColor());
                drawPoint(g2, p);
            }
        }

        public void drawLines(Graphics2D g2)
        {
            clearFlags(HalfEdge.FLAG_DRAW);
            foreach (HalfEdge he in halfEdges)
            {
                g2.setColor(he.getColor());
                if (DRAW_HALFEDGES) drawHalfEdge(g2, he);
                if (he.isFlagged(HalfEdge.FLAG_DRAW)) continue;
                drawLine(g2, he);
                he.flagEdge(HalfEdge.FLAG_DRAW);
            }
        }

        public void drawLinesGray(Graphics2D g2)
        {
            clearFlags(HalfEdge.FLAG_DRAW);
            g2.setColor(grayColor);
            foreach (HalfEdge he in halfEdges)
            {
                if (he.isFlagged(HalfEdge.FLAG_DRAW)) continue;
                drawLine(g2, he.origin, he.next.origin);
                he.flagEdge(HalfEdge.FLAG_DRAW);
            }
        }

        public void drawBoundary(Graphics2D g2)
        {
            g2.setColor(HalfEdge.COLORS[(HalfEdge.BOUNDARY)]);
            foreach (HalfEdge he in halfEdges)
            {
                if (he.isType(HalfEdge.BOUNDARY)) drawLine(g2, he);
            }
        }

        protected void drawPoint(Graphics2D g2, Point2d p)
        {
            if (DRAW_PT_FULL)
            {
                g2.fillOval((int)p.X - 3, (int)p.Y - 3, 5, 5);
                g2.drawOval((int)p.X - 6, (int)p.Y - 5, 10, 10);
            }
            else
            {
                g2.fillOval((int)p.X, (int)p.Y, 1, 1);
            }
            if (DRAW_PT_INDICES)
            {
                g2.drawString(
                        String.valueOf(points.IndexOf(p)),
                        (int)p.X + 5,
                        (int)p.Y + 5
                        );
            }
        }

        protected void drawLine(Graphics2D g2, HalfEdge he)
        {
            Point2d p1 = he.origin;
            Point2d p2 = he.next.origin;
            g2.drawLine((int)p1.X, (int)p1.Y, (int)p2.X, (int)p2.Y);
            if (DRAW_HE_INDICES)
            {
                int x1, y1, x2, y2;
                double x, y;
                Vector2d v = new Vector2d(p2.X - p1.X, p2.Y - p1.Y);
                v.normalize();
                v.scale(5);
                x = v.X - v.Y;
                y = v.X + v.Y;
                x1 = (int)(p1.X + x);
                y1 = (int)(p1.Y + y);
                x2 = (int)(p2.X + x);
                y2 = (int)(p2.Y + y);
                g2.drawString(String.valueOf(halfEdges.IndexOf(he)),
                        (int)(0.5 * (x1 + x2)),
                        (int)(0.5 * (y1 + y2)));
            }
        }

        protected void drawLine(Graphics2D g2, Point2d p1, Point2d p2)
        {
            g2.drawLine((int)p1.X, (int)p1.Y, (int)p2.X, (int)p2.Y);
        }

        protected void drawHalfEdge(Graphics2D g2, HalfEdge he)
        {
            drawHalfEdge(g2, he.next.origin, he.origin);
        }

        protected void drawHalfEdge(Graphics2D g2, Point2d p1, Point2d p2)
        {
            int x1, y1, x2, y2;
            Vector2d v = new Vector2d(p2.X - p1.X, p2.Y - p1.Y);
            v.scale(0.1);
            x1 = x2 = (int)(p1.X + v.X);
            y1 = y2 = (int)(p1.Y + v.Y);
            v.normalize();
            v.scale(4);
            x2 += (int)(v.X + v.Y);
            y2 += (int)(v.Y - v.X);
            g2.drawLine(x1, y1, x2, y2);
        }
         */
        #endregion

        /*
        protected string error(String s, Object ... args)
        {
            if (DEBUG) {
                debugPause(s,args);
            }
            return error(s);
        }
         */
    }
}